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Abstract 

For the calculation of neutral excitations, time-dependent density functional theory (TDDFT) 
is an exact reformulation of the many-body time-dependent Schrodinger equation, based 
on knowledge of the density instead of the many-body wavef unction. The density can 
be determined in an efficient scheme by solving one-particle non-interacting Schrodinger 
equations — the Kohn-Sham equations. The complication of the problem is hidden in the — 
unknown — time-dependent exchange and correlation potential that appears in the Kohn-Sham 
equations and for which it is essential to find good approximations. Many approximations have 
been suggested and tested for finite systems, where even the very simple adiabatic local-density 
approximation (ALDA) has often proved to be successful. In the case of solids, ALDA fails to 
reproduce optical absorption spectra, which are instead well described by solving the Bethe- 
Salpeter equation of many -body perturbation theory (MBPT). On the other hand, ALDA can 
lead to excellent results for loss functions (at vanishing and finite momentum transfer). In 
view of this and thanks to recent successful developments of improved linear-response kernels 
derived from MBPT, TDDFT is today considered a promising alternative to MBPT for the 
calculation of electronic spectra, even for solids. After reviewing the fundamentals of TDDFT 
within linear response, we discuss different approaches and a variety of applications to extended 
systems. 

(Some figures in this article are in colour only in the electronic version) 
This article was invited by professor M W Finnis. 
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1. Introduction 



Most of today's quantum mechanical theoretical research in condensed matter physics and 
chemistry is not aimed at finding new fundamental interactions or basic laws — it deals 
with solving the Schrodinger equation of a well-known Hamiltonian and extracting useful 
information from the solution. This Hamiltonian, however, describes a many-body problem, 
and for a number of electrons well above 10 it is impossible to even dream of a full numerically 
exact solution. Moreover, the exact solution would yield a wealth of information that could 
hardly be understood without further analysis and simplification and would contain many 
details that, for a given situation or question, one is probably not interested in [1]. Therefore, 
it is often more appropriate to reformulate the problem from the start, working with effective 
Hamiltonians or selected expectation values that are suitable for the solution of a reduced 
problem. This procedure will ideally simplify both the calculation and the analysis of the 
desired quantities. 

Density-functional theory (DFT) [2] is a prominent example of such an approach. It has 
been designed for the calculation of ground- state properties. It is based on knowledge of the 

density n(r) instead of the full many -body wavefunction *I>(ri, r 2 , , r N , o\, a 2 , cr N ) of 

the N -particle system. DFT can be formulated in the Kohn-Sham (KS) approach [3] where 
an efficient one-particle non-interacting Schrodinger equation — the Kohn-Sham equation — 
yields eigenvalues s t and orbitals <Pi(r). The orbitals and eigenvalues do not in general have 
a direct physical meaning, but the former can be used to construct the true density of the 
interacting system according to n(r) = J^. |<p;(r)| 2 . The complication of the problem is 
hidden in the — unknown — exchange and correlation (xc) total energy E xc and the exchange- 
correlation potential v xc [n](r) that appears in the Kohn-Sham equation. Very efficient 
approximations have been proposed, such as the local-density approximation (LDA) [3] or 
generalized gradient approximations (GGA) [4,5], and many ground-state properties are today 
calculated from first principles with a precision of a few per cent, such as lattice parameters 
or phonon frequencies [6]. There exist, however, ground- state properties for which even 
in simple systems standard approximations do less well: cohesive energies in particular can 
easily be off by 10% in LDA (because of errors in calculating the isolated atoms that enter this 
total-energy difference), and failures are reported for static response properties, such as the 
dielectric constant €oo, which is often substantially overestimated [7]. Other problems arise, 
e.g. in the description of strongly correlated systems [8] or of the van der Waals dispersion 
attraction [9] . These problems in calculating ground-state properties can be traced back to the 
limits of validity of the employed approximations. 

Another problem of static ground-state DFT-Kohn-Sham is the fact that excitations, such 
as those measured in the optical response to a time-dependent electric field, are in principle 
not accessible. This is not a question of the available approximations but of the fact that the 
theory is not meant to describe these phenomena. In fact, even if one could calculate the exact 
Kohn-Sham eigenvalues, their differences would not necessarily be close to the measured 
excitation energies. By definition, they do not stand for electron addition or removal energies 
either [10] . The fact that the Kohn-Sham gap is in general reported to be too small with respect 
to the measured gaps therefore does not a priori tell us anything about the quality of a chosen 
approximation for the exchange-correlation potential. 

If one wishes to work with an efficient Hamiltonian that in principle yields eigenvalues 
meant to be electron addition or removal energies, or excitation energies, more than just the 
static ground-state density has to be calculated. Such energies can be found essentially in 
two ways. 
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First, by studying particle propagation and fluctuations in the system. This yields 
correlation functions that can then be related to response functions yielding, e.g. linear response 
for optical absorption. These correlation functions are one- or two-body Green's functions [11] 
(or higher-order ones, for problems beyond the scope of this review). The one-body Green's 
function (that can essentially be understood as a time-dependent particle and hole density 
matrix) has phase fluctuations (or, in frequency (co) space, poles) given by electron addition 
and removal energies measured, e.g. in photoemission or inverse photoemission experiments. 
The particle-hole part of the two-particle Green's function, in turn, has poles at the energies 
of neutral excitations. A contraction of the four-point reducible part L(r,r\,r f ,r[, co) of 
the two-particle Green's function leads to the two-point response functions x( r > r > <*>) mat 
determine measurable spectra, such as absorption or electron energy loss spectra (EELS). 
Many-body perturbation theory (MBPT) yields a framework where suitable approximations for 
these Green's functions can be found; in particular, the Bethe-Salpeter equation (BSE) is a good 
starting point for approximations for x [11-14]. The price to be paid for a physically intuitive 
and in general quite reliable description is however relatively high in terms of computational 
cost, because now quantities such as L(r, r\,r\ r[, co) appear, instead of the density n(r). 

Second, by actually (on paper or using a computer) exposing the system to a time- 
dependent external potential and calculating the evolution of the density in time. The response 
function x, for example, can then be directly determined from the linear-response relation 
ft (1) (r, co) = f dVx(f% r r , co)v^ t (r / , co) between the variation in the external potential and 
the induced density. This route has become accessible owing to the extension of DFT to its 
time-dependent generalization, TDDFT [15-19]. Put on a rigorous basis by the Runge-Gross 
theorem [16], one can understand that in TDDFT the quantum-mechanical 'trajectory' of the 
system under the influence of a time-dependent external potential is found by searching for the 
extrema of an action (instead of the minimization of a total energy, as done for the ground state), 
by analogy to the case of classical mechanics. One obtains therefore the time-dependent Kohn- 
Sham equations as a generalization of the static case and from these the response functions 
describing neutral excitations of a system [20] . At this point the difficulty resides in finding 
suitable approximations for the time-dependent exchange-correlation potential v xc [n](r, t)\ 
note that now the functional dependence is on the density in the whole space and at all past 
times. 

Many approximations have been suggested and tested for finite systems. Even the 
very simple adiabatic local-density approximation (ALDA, also called time-dependent LDA 
(TDLDA)) where v^ DA [n](r, t) = v^ A (n(r, t)) has been proved to be very successful in 
many cases [15, 21], although the lack of a long-range (1/r) decay of the potential can lead to 
serious problems for questions such as Rydberg states [22]. The latter shortcoming of the LDA 
potential is not so crucial in solids where the electron density is quite homogeneous (compared 
with an atom in empty space); instead, the wrong long-range behaviour of the linear-response 
kernel f xc (T,r' t, t') = 8v xc [n](r, t)/8n(r f , t') can cause serious problems [14]. In fact, in 
the ALDA this kernel is proportional to 8(r — r f ), whereas in non-metallic systems it should 
decay as l/\r — r f \ [23]. This shortcoming already shows up, for example, in the calculation 
of polarizabilities for molecular chains [24]. In the case of absorption spectra of solids, where 
the imaginary part of the dielectric function € for the vanishing wavevector q (corresponding 
to a macroscopic average) is calculated, the lack of a l/q 2 divergence (stemming from the 
Fourier transform of l/|r — r'\) can lead to drastic failures. For example, the ALDA is not 
able to reproduce bound excitons [14]. On the other hand, for finite momentum transfer or 
when the loss function — Im(€ -1 ) is the quantity of interest (e.g. in electron energy loss or 
inelastic x-ray scattering spectra) this term is not dominant, and ALDA can lead to very good 
results (see, e.g. [25-28]). For this reason, and because of recent successful developments of 
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improved linear-response kernels derived from MBPT [29-35], TDDFT is today considered 
to be a promising alternative to MBPT for the calculation of electronic spectra, even for solids. 

This situation is the main motivation for summarizing here the state-of-the art of the 
TDFFT approach for extended systems. This review is centred on the calculation of the linear 
response, which actually constitutes the overwhelming majority of work done in the field. 
Moreover, it is limited to the calculation of excitation of valence electrons. TDDFT is also used 
to determine core electron spectra (see, e.g. [36]); however, this case is in many respects more 
similar to the situation in finite systems and will therefore not be treated here. Complementary 
to the present review, a discussion of time-dependent current density functional theory for 
periodic systems can be found in [37]. A comprehensive review of many aspects of TDDFT 
can be found in [38]. 

After a brief review of the fundamentals of TDFFT for which several more detailed reviews 
can be consulted (see, e.g. [19,21,37]), we will concentrate on questions that are more specific 
to extended systems, with a short comparison with finite systems. Rather than giving an 
exhaustive summary of all results that have been obtained in this rapidly expanding field, 
selected applications are presented to illustrate major points, such as the importance of crystal 
local-field effects or the effect of contributions derived from MBPT. A short outlook follows 
the conclusions. 

2. Fundamentals 

As the quantum-mechanical treatment of stationary and time-dependent systems differs in 
many aspects, it is not straightforward to generalize the mathematical framework of static 
density-functional theory. For example, the total energy, which plays a central role in the 
original Hohenberg-Kohn theorem [2], is not a conserved quantity in the presence of time- 
dependent external fields, and there is hence no variational principle for it on the basis of the 
density that can be exploited. In this section we start by discussing the theoretical foundations 
of TDDFT with a special emphasis on the linear density -response function and its connection 
to the electronic excitation spectrum. 

We will use atomic units throughout the paper (i.e. e 2 = h = m e = 1). 

2.1. The Runge-Gross theorem 

The evolution of a (non-relativistic) spin-unpolarized interacting many-electron system is 
governed by the time-dependent Schrodinger equation 



where H is the Hamiltonian operator of the system and {r} = {r\, . . . , r N } are the spatial 
coordinates of the TV electrons. The Hamiltonian can be written in the form 



where v ext (r, t) is the time-dependent external potential and vfo — rj) = l/|r; — rj\ the 
Coulomb interaction. Being interested in spectroscopy, we consider scenarios where the 
system is initially at rest in a static potential v ext (r, t) = v® xt (r), before a time-dependent 
perturbation is switched on at t = to in order to probe the response of the electron system. Under 
these circumstances the initial state at to is given by the stationary ground- state wavef unction 
*o) = ^o(W) exp(— iEoto), where Eq denotes the ground-state energy. This initial 
wavef unction is determined up to an irrelevant phase factor for non-degenerate systems. 



i-V({r},t) = H({r},t)V({r},t), 



V({r},t ) given, 



(2.1) 




(2.2) 
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By virtue of the Hohenberg-Kohn theorem it is also a functional of the static ground-state 
density n(r, to) = hqs (**). We only admit physical potentials that are finite everywhere and 
vary smoothly in time, so that they can be expanded into a Taylor series about the initial time to'. 



t(r ' ^ = J2 Cj rr^ ~ ^ k with Ckir) = ^ ext(r ' f) 

k=0 



(2.3) 



The theoretical basis of TDDFT is the Runge-Gross theorem [16], which asserts the one-to- 
one correspondence between the external potential and the density, thus playing the same role 
as the Hohenberg-Kohn theorem in static density-functional theory. Of course, for a given 
external potential it is always possible, in principle, to solve the time-dependent Schrodinger 
equation (2.1); the density is then given by 

n(rJ) = N jd 3 r 2 jd 3 r 3 • • • jd 3 r N \V(r, r 2 , . . . , r N , t)\ 2 . (2.4) 

What remains to be proved, in order to demonstrate the one-to-one correspondence, is that if 
two potentials v ext (r, t) and v f ext (r, t) differ by more than a purely time-dependent function, 
then the associated densities n(r, t) and n\r, t) must be distinct. The addition of a purely 
time-dependent function is exempt because it only changes the phase of the wavefunction 
but not the density. One assumes that both systems evolve from the same initial ground-state 
wavefunction ^ ({r} , to) . The expansion coefficients of the two potentials around to are denoted 
by Ckir) and c r k (r), and one defines Uk(r) = Ck(r) — c' k (r). If the potentials differ by more than 
a purely time-dependent function, then at least one coefficient Uk(r) is not a mere constant 
but a spatially varying function. For the proof of the Runge-Gross theorem given in [16] one 
makes use of the current density 

j(r, t) = - l -N j d 3 r 2 J d 3 r 3 - J d 3 r N {^*(r, r 2 , . . . , r N , f)[V*(r, r 2 , . . . , r N , t)] 

-[V^/*(r, r 2 , . . . , r N , t)]V(r, r 2 , . . . , r N , t)} , (2.5) 
which can also be written in a second quantization formalism as 

j(r, t) = -i<*(f)l#V)[V#(r)] - [V#V)]#(r)|*(f)). ( 2 - 6 ) 
The time evolution of the current density can be discussed by means of the equation of motion 

i^/(r, t) = (*(f)l[/V), ff(f)]|*(f)). (2.7) 
at 

Moreover, j(r, t) is related to the density through the continuity equation 

%-n(r,t) = -V-j(r,t). (2.8) 

Ot 

This identity expresses the conservation of the total particle number in a differential form: the 
change in the number of electrons within a certain volume equals the flux through its surface. 
In the first step one shows that the current densities j(r, t) and f(r, t) induced by the two 
potentials differ. To this effect one examines the time derivative 

i^{/(r, t) -f(r, t)} t=t0 = (V \m, H(t ) - ff'fo)]|*o> 

= (^o\\j(r), fiatfr, t ) - C; xt (r, t )]\Vo) 

= in(r, t )V{v ext (r, t ) - v' ext (r, t )} = n(r, t )Vu (r), (2.9) 
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which follows from the definition (2.6) together with the known evolution of the current density 
(2.7). If uo(r) is not a constant, then the right-hand side is non-zero, and consequently the 
derivatives of the current densities at to must be distinct. The potentials might also differ by a 
coefficient Uk with k ^ 0; in this case one can take an appropriate higher time derivative 

^{/(r, t) -fir, t)} t=t0 = n(r, fo)VKjfc(r), (2.10) 

with the non-constant u k ir) hence establishing that at least one term in the Taylor expansions of 
j(r, t) and f(r, t) differs. This implies that the current densities themselves deviate for t > to. 
In the second step one proves that the corresponding densities also differ. For this purpose 
one takes the (&+l)st time derivative of the continuity equation (2.8) and again examines the 
difference 

d k + 2 d k + l 

HkzW* f) ~ n <?> = " v ' sp*iV<r> t) -fir, t)} t=t0 

= -V.{nir,t )Vu k ir)} (2.11) 

between the two systems. If the quantity on the right-hand side is non-zero, then the (£+2)nd 
terms in the Taylor expansions of n(r, t) and n'ir, t) around to differ, and the densities 
themselves must hence deviate for t > to. The original proof [16] refers only to finite systems, 
where both the potential and the density decay to zero at large distances, but for extended 
systems it is easy to see that the right-hand side of (2.1 1) vanishes only if Uk(r) is of the form 

C r n(0, to) 

u k ir) = u k i0)- -±^-E k -dr / , (2.12) 
Jo nir'Jo) 

with an arbitrary but constant vector E k . As the density is always positive, u k ir) then grows 
beyond all bounds as \r\ -> oo, which implies that at least one of the potentials v ext (r, t) or 
v f ext (r, t) becomes infinite. However, this case was explicitly excluded. For all finite physical 
potentials the right-hand side of (2.11) is indeed non-zero. This concludes the proof of the 
Runge-Gross theorem. Note that potentials of the type (2.12) are also incompatible with the 
Hohenberg-Kohn theorem in static density-functional theory: as the energy of the electrons can 
always be lowered by a translation in the direction of the field vector, there is no ground- state 
solution [39]. 

The Runge-Gross theorem is, in fact, a very strong statement: from knowledge of 
the density alone it is possible to deduce the external potential and hence the many-body 
wavefunction, which in turn determines every observable of the system. Therefore, all 
observables can ultimately be regarded as functionals of the density. We note that, in contrast 
to more general cases [40,41], there is no additional initial-state dependence in this scenario, 
because the stationary wavefunction at to itself is determined by the static ground- state density 
?7 GS (r) = nir, to) of the unperturbed system. 



2.2. The time -dependent Kohn-Sham equations 

The Runge-Gross theorem states that all observables are functionals of the density, but it 
contains no prescription on how this central quantity can actually be calculated. To overcome 
the analogous problem in static density-functional theory, Kohn and Sham [3] suggested to 
use an auxiliary system of non-interacting electrons moving in an effective local potential, 
which is designed in such a way that the densities of the non-interacting system and the 
real interacting electrons coincide. This scheme has the big advantage to include the exact 
non-interacting kinetic energy, which represents almost all the true kinetic energy of the N- 
electron system. The main task is then to find a good approximation for this a priori unknown 
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effective potential. This idea was generalized to the time-dependent case, where the Kohn- 
Sham electrons obey [16] 

i^<Pj(T, t) = (~V 2 + UKsW(r, t)j <pj(r, t) (2.13) 

and the density is given by 

N 

n(r,t) = J2\<Pj(r>t)\ 2 - ( 2 - 14 ) 

7 = 1 

The Kohn-Sham scheme assumes that one can always find a local potential v^[n](r,t) 
with the property that the orbitals obtained from (2.13) reproduce the given density of an 
interacting electron system, but the validity of this assumption, known as 'non-interacting 
u-representability', is not obvious and requires a careful examination. If such a potential 
exists, however, then by virtue of the Runge-Gross theorem it is unique up to a purely time- 
dependent function. Giving a constructive proof, van Leeuwen [42] showed that an effective 
local potential with the desired property exists if one can find a stationary wavefunction that 
yields the initial density n(r, to) and is the ground state of a non-interacting electron system. 
The problem is thus reduced to the question of non-interacting u-representability in static 
density-functional theory. Despite much progress, the latter is still unresolved. Examples of 
well-behaved densities that do not correspond to the ground state of a non-interacting system 
are known [43,44]; the implications of this discovery remain however unclear. In actual 
calculations, where the initial Kohn-Sham wavefunction is obtained from the constrained 
minimization of a smooth approximate energy functional, a solution can always be found [45]. 

If non-interacting f-representability is assumed, then VKsW(r,t) is determined 
completely by the requirement that (2.14) equals the density of the real interacting electron 
system. As in the case of ground-state DFT, one has then to find an explicit expression for the 
effective potential that can be exploited to construct useful approximations. For this purpose 
it is convenient to employ the same separation 

VKsW(r, t) = v QXt (r, t) + v n [n](r, t) + v xc [n](r, t) (2.15) 

as in static density-functional theory. The first term is the external potential, the second is the 
Hartree potential 

v H [n](r,t)= /dV-^ (2.16) 



and the third incorporates all remaining exchange and correlation effects. In the static case one 
can exploit the variational principle and determine the orbitals of the Kohn-Sham electrons in 
such a way that the total energy is minimized; all potential terms are then obtained as functional 
derivatives of the corresponding energy contributions with respect to the density. The energy 
in turn is a well-defined physical quantity and amenable to approximations. In systems driven 
by time-dependent external fields the total energy is not a conserved quantity and there cannot 
be a minimization principle. There exists, however, a quantity analogous to the energy, the 
quantum-mechanical action functional 

AW = jf 1 dt(d>(t)\ (iy t - H(f)j |<D(0> , (2.17) 

which has the property that its derivative with respect to a Af-body function (O (t) | vanishes at 
the true many-body wavefunction, i.e. the solution of the Schrodinger equation 

8Am =[ i T<- H(t) ) !* w > = °- (2 ' 18) 



5(O(0l 



Time-dependent density-functional theory for extended systems 



365 



C 



t t h 

Figure 1. The Keldysh contour C, starting at to and turning back at t\ . Pseudotime values t on the 
forward and backward branches are distinct. 



Therefore, it is possible to solve the time-dependent problem by searching for the stationary 
point of the action. In contrast to the energy in the static case, the stationary point is not 
necessarily a minimum, however. Furthermore, the value of the action itself does not provide 
any relevant additional information, since for the true many -body wavefunction A[^] =0. 

By virtue of the Runge-Gross theorem we may consider the action as a functional of the 
density. The obvious definition of A[n] is to evaluate (2.17) at the wavefunction 0[ft]({r}, t) 
that evolves from the given initial state and yields the density n(r,t). In analogy to the total 
energy in static density-functional theory, one would expect that the true density makes this 
functional stationary and can thus be identified. A suitable decomposition of the action would 
then define the exchange-correlation potential in terms of the functional derivative 
8A xc [n] 

Vn(r,t) = —p-±. (2.19) 
8n(r, t) 

Unfortunately, this procedure is doomed to failure. A first problem arises because the 
density determines the potential only up to a purely time-dependent function. Therefore, 
the wavefunction 0[/t]({r}, t) and the value of the action derived from it are not unique. Even 
if the phase of the wavefunction is fixed by imposing additional constraints, there is another 
more fundamental problem, which becomes evident if one examines the second functional 
derivative 

8vxc(r, t) _ 8 2 A xc [n] 

8n(r f ,t') ~ 8n{r,t)8n(r',t')' 
Whereas the expression on the right-hand side is symmetric in (r, t) and (r f , t f ), the exchange- 
correlation potential can only be influenced by the density at earlier times. Therefore, causality 
dictates that the left-hand side must vanish for t < t' but not for t > t' . The symmetry and 
causality requirements contradict each other and cannot be satisfied simultaneously. One is 
hence forced to conclude that a differentiable functional A xc [n] with the property (2.19) does 
not exist. 

This causality dilemma [46, 47] was eventually resolved by van Leeuwen [45] using the 
time-contour method given by Keldysh [48] . In this approach to non-equilibrium dynamics the 
physical time t is parametrized by an underlying parameter r called pseudotime in such a way 
that t(x) runs from to to t\ and back to to if r runs along the contour C illustrated in figure 1. 
As pseudotime values on the forward and backward branches are distinct, the ordering along 
the contour differs from that on the physical time axis. The solution of the dilemma hence 
consists of satisfying the causality and symmetry requirements in different variable spaces. To 
this effect the action is first defined as a functional of the external potential in a form that does 
not explicitly contain 3/3 1: 

A[U]=i\n(V(to)\T c exp(-i f drt , (r)H(r)\ |*(f )> . (2.21) 



For the derivation of (2.21) the reader can see [45]. The potential U(r,r) is contained in 
the Hamiltonian, and Tq sorts the subsequent operators in the order of ascending pseudotime 
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arguments from right to left. For physical potentials of the form U(r,z) = v ext (r, t(z)) the 
value of the action is zero, because the contributions along the two time-contour directions 
cancel each other, but its derivative can be non-zero. In fact, the action (2.21) is defined in 
such a way that its functional derivative yields the density 

8A[U] 

= n(r,z). (2.22) 



8U(r, z) 

An unambiguous functional of the density can then be constructed by means of a Legendre 
transformation 

A[n] = -A[U] + j dzt'(z) j d 3 rU(r, z)n(r, z). (2.23) 

Finally, for practical purposes the action is decomposed according to 

If , f o f o f n(r, z)n(r f , z) 
A[n] = A KS [n] --J dzt'(r) J d 3 r J dV ' '_ K ' - A xc [n]. (2.24) 

The first term is the action of the non-interacting Kohn-Sham system, whose Legendre 
transform A K s[^ks] is defined in analogy to (2.21) in terms of the initial Kohn-Sham 
wavefunction and the effective local potential. The second term is related to the Hartree 
potential and the third gives rise to the exchange-correlation potential 

8A xc [n] 



Vxc(r, t) 



8n(r, z) 



(2.25) 

n(r,t) 



Defining the action in the pseudotime domain instead of the real time axis guarantees the 
proper symmetry of the second functional derivative in (r, z) and (r r , z'). On the other hand, 
the exchange-correlation potential (2.25), which is obtained by inserting the physical time 
argument after performing the functional derivative with respect to n(r, r), respects causality 
on the time axis. From a theoretical point of view, all quantities that enter the Kohn-Sham 
scheme are thus well defined, and working approximations can be derived by finding a suitable 
expression for the action functional, for example, through an expansion in the powers of the 
Coulomb interaction. This approach is known as the time-dependent optimized effective- 
potential method [49]. Unfortunately, the leading term, which is linear in the Coulomb 
interaction and retains only exchange and no correlation [50], has already a high computational 
cost. In fact, at present the design of specific approximations for the exchange-correlation 
potential in TDDFT is still at a very early stage; we will discuss some promising approaches 
later in this review. Today, most calculations take however a pragmatic point of view and 
simply use one of the established functionals of static density-functional theory. The most 
popular choice is the adiabatic local-density approximation (ALDA) [15], which is obtained 
by evaluating the standard LDA potential with the time-dependent density n(r,t): 

v£° A [n](r, t) = v™ G [n](r, t)\ n=nM . (2.26) 

The adiabatic approach is a drastic simplification, however, and a priori only justified for 
systems with a weak time-dependence which are always locally close to equilibrium. This 
adds to the problems that are due to the spatial locality of the LDA. 



2.3. Linear-response theory 

If the time-dependent external perturbation in v ext (r, t) = v^ t (r) + v^ t (r, t) is weak, then 
linear-response theory can be exploited to describe the dynamics of a system more efficiently 
than a full solution of the Kohn-Sham equations (2.13). In this case the density is expanded 
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in orders of v^ t (r, t) according to n(r, t) = n^(r) + n^ l \r,t) + • •, where the first-order 
correction is given by 

n (1) (r,0 = j dt f jd 3 r' X (r,r\t (2.27) 

in terms of the linear density-response function 

y(r,r ,t — t ) = 

Causality requires x( r ^ r ^ ~ t') = for t < f, of course, because the density cannot be 
influenced by later variations of the potential. To calculate the linear density -response function 
in practice one exploits the fact that the density of the real system is identical to that of the 
non-interacting Kohn-Sham electrons. As the latter move in the effective potential v^sir", t"), 
one starts by applying the chain rule for functional derivatives 



(2.28) 



x (r, /,*-*') = I™ At" fdV 



8n(r,t) 8v KS (r",n 
8vKs(r»,t») 8v ext (r',t>) ' 



(2.29) 



The first term on the right-hand side corresponds to the linear density-response function 
Xks(t% r" ,t — t") of the non-interacting Kohn-Sham system, since the effective potential plays 
the role of the 'external potential' of the Kohn-Sham system. It can be calculated explicitly 
from time-dependent perturbation theory and is given by 

XKs(f, r , co) = hm > > (fj - /*)— ; (2.30) 

1^0+ ff— f CO- S k +Sj +IY) 

7 = 1 k=\ J 

in frequency space. The energies Sj appearing in the denominator are the eigenvalues 

corresponding to the unperturbed stationary Kohn-Sham wavefunctions cpjir). In order to 

evaluate the second term in (2.29) one uses the separation (2.15), which yields 

8v KS (r",t") , „ , 8v ¥i (r"j") 8v xc (r\t") 

SV } = 8{r" - r')8(t" - t') + V ' + cV - . (2.31) 
$v Q Ar',n 8v QXt (r',t>) 8v QXt (r',t>) 

As both the Hartree potential and the exchange-correlation potential are functionals of the 

density, one can apply the chain rule once more and rewrite these two contributions as 



$VcAr f ,t>) J r 8n(r'",t">)8v QXi (r>,t>) 



and analogously for 8v xc (r", t") / 8v QXt (r f , t'). The last term on the right-hand side of (2.32) 
is easily recognized as the linear density-response function x( r " f i r > t'" — f). The functional 
derivative of the Hartree potential with respect to the density follows from the definition (2.16) 
and simply equals the Coulomb potential 

Sv »™=-!— S(t »-n. (233) 
8n(r"',t"') \r" -r"'\ 

The last term of (2.31) contains the so-called exchange-correlation kernel 

8v xc (r",t") 



f xz ir",r'\t" -*"') 



(2.34) 

n(r"' \t"')=n {0 \r"') 



8n(r f ", t>") 

After collecting all terms and performing a Fourier transform to frequency space, whereby 
convolutions on the time axis turn into simple multiplications, one obtains the final integral 
equation [20] 

X(r, r\ co) = X Ks(r, r\ co) + j d 3 r" j dV" X Ks(r, r\ co) 

x + U c {r\r>\co)^x{r"\r\co). (2.35) 
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The TDDFT equations in the linear-response regime can be cast in numerous different 
forms. For solids, in most implementations the integral equation (2.35) is solved routinely by 
projecting all quantities onto a suitable set of basis functions. Very often, one uses a plane wave 
representation within the pseudotential approximation (see e.g. [51,52]), but localized basis 
sets can equally be used to allow for all-electron calculations (see e.g. [53]). Equation (2.35) 
thus turns into a matrix equation /(&>) = Xks(&0 + Xks(<^)[^ + /xc(&0]x(<^)> f° r example, in 
reciprocal lattice vectors in the case of periodic systems where x = Xg,g '((Z)- If one wishes 
to only obtain an absorption spectrum or the loss function at a given momentum transfer (see 
definitions in section 3) only one component of the matrix xg g <*>) is required. This can be 
obtained by solving a linear system of equations, thus avoiding a numerically involved matrix 
inversion [51]. 

Alternatively, absorption spectra can be calculated by propagating the full TD Kohn- 
Sham equations in real time [54]. This description decreases storage requirements; it allows 
the entire frequency-dependent dielectric function to be calculated at once, and the scaling 
with the number of atoms is quite favourable. However, the prefactor is fairly large as such 
calculations typically require ^10000 time-steps with a time-step of ~10 -3 fs [55]. 

Another efficient approach, based on the linear response within Ghosh and Dhara's time- 
dependent density-functional formalism [56], was proposed [57]. It uses an iterative scheme 
in real space, in which the density and the potential are updated in each cycle, thereby avoiding 
the explicit evaluation of the Kohn-Sham response kernels. 

Finally, a method to calculate the dynamical polarizability using only occupied states 
has been proposed recently [58]. The dynamical polarizability is represented by a matrix 
continued fraction whose coefficients can be obtained from a Lanczos method. This method 
scales favourably with the system size, and it may become useful for large scale systems. At 
present there is, however, only a single application to the benzene molecule [58]. 

2.4. Excitation energies 

In static DFT the interpretation of the one-particle Kohn-Sham eigenvalues Sj as quasiparticle 
energies is not formally justified and it leads to the well-known problem of the underestimation 
of transition energies. In the framework of TDDFT the relevant information about the excited 
states is contained in the linear density-response function: in fact, it can be shown that the 
true excitation energies are the poles of x 0% r > <*>) • In contrast to other attempts to calculate 
electronic excitations within a density-functional framework [59-61], TDDFT has the great 
advantage that it is not restricted to a subset of excited states but, in principle, yields the 
complete excitation spectrum. 

In order to see this, one can calculate the density change due to the external potential at 
first order. The stationary eigenstates of the original unperturbed Hamiltonian are labelled by 
tyj({r], t) = ^ ({r}) exp(— iEjt), where Ej denotes the corresponding energy eigenvalues. 
After the onset of the time-dependent perturbation it is possible to expand the wavefunction 
^({r}, t) that evolves from the ground-state ^ (0) ({r}, t) = ^(M) exp (— i£oO in orders of 
v^ t (r, t). The first-order correction is 





(2.36) 
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The corresponding change in the density is 

.<»<,.,) -*/dV ./<.%([*<»<,,, r».or«»fr.r, r».0 

+[^ (0) (r, r 2 , . . . , r#, r 2 , . . . , r N , t)}. (2.37) 

In order to simplify the notation we introduce the overlap functions 

rijir) = N j d 3 r 2 • • • J d 3 r N V*(r, r 2 , . . . , r N )^j(r, r 2 , . . . , r N ), (2.38) 

and after inserting (2.36) and (2.38) into (2.37) we obtain 



J —oo J 



n^(r,t) = df d 



-i^(n y -(r)n*(/)e- i( ^ 



-E )(t-t') 



- nJ(r)/i ; -(Oe i(£r£o)(f "°)0(r - *')] (2.39) 

Comparing this expression with (2.27), one finds that the term in square brackets equals the 
linear density-response function x (r, ^, t — t r ). The Heaviside step function ®(t — t f ) has been 
introduced to replace the integral over time dt' with f™ dt f . After a Fourier transform to 
frequency space and using ®(t) = i/27r lim^o+ ^ C0 '^j Q ~ lta) one arrives at the Lehmann 
representation of the density -response function: 

X(r, r', co) = lim V ( 7 J J 7 ) , (2.40) 

^^o+ ^— ^ \co — E j + Eq + IT] CO + E j — Eq + IT] J 

where ?7 is a positive infinitesimal. From (2.40) it is evident that the poles of x( r > r > &>) 
correspond to the exact excitation energies Ej — Eq. Furthermore, all quantities on the right- 
hand side depend only on the Hamiltonian of the unperturbed stationary system. By virtue of 
the Hohenberg-Kohn theorem the linear density-response function is hence a functional of the 
static ground- state density ncs( r )- 

The form of (2.40) is valid for finite systems with discrete eigenvalues. As the energies Ej 
of the eigenstates of the many-electron system are real, it appears that the poles of x (r, r' , co) 
are at real energies. For extended systems, on the other hand, the spectrum is continuous, and 
the sum in (2.40) turns into an integral that gives rise to a branch cut along the real energy axis. 
The infinitely close-lying resonances thus merge into broad structures that can be identified 
with elementary quasiparticles, such as plasmons or excitons. As these structures have a certain 
width, they are described by poles in the complex plane with a real part, which corresponds to 
the energy of the excitation, and a finite imaginary part, whose inverse is proportional to the 
excitation lifetime. 



3. TDDFT in practice: approximations and problems 

Linear-response theory can be applied now to study the response of an extended system to a 
small time-dependent perturbation v QXt (r, t). The linear variation in the density induced by the 
perturbation is given by (2.27). As a consequence of the polarization of the system due to the 
applied perturbation, the total potential becomes a sum of the external potential and the induced 
potential: v tot = v QXt + v^. The basic quantity that gives information about the screening of 
the system in linear response is the microscopic dielectric function €, which relates the total 
potential v tot to the applied potential v QXt : 

^(r, t) = [ dt' [ dVc-Vy, t - t')v ext (r\ f). (3.1) 

J —oo J 
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The microscopic dielectric function € and the reducible polarizability x are hence related by 

€~ l (r,r\ t -t') = 8(r-r f )8(t - t') + j dV 'v(r - r")x(r" \r> \ t - t'). (3.2) 

For periodic systems, the most natural way to deal with spatial periodicity is to apply a Fourier 
transform and rewrite (3.2) in reciprocal space: 

€~^(q, co) = 8gg> + v G (q)xGG'(q, o>), (3.3) 

where G is a vector of the reciprocal lattice, while q is a vector in the first Brillouin zone. In 
(3.3) a Fourier transform has also been applied to move from time to frequency space. 

From the microscopic dielectric function one has to obtain measurable quantities. In 
the case of absorption spectra, this means to calculate the imaginary part of the macroscopic 
dielectric function [62-64] : 

1 

€ M (co) = hm - — =— . (3.4) 

<7-o [€ G l G ,(q, co)] GiG > =0 

We have dealt so far with the response to a potential whose electric field is longitudinal with 
respect to the wavevector. Light, instead, is a transverse perturbation, i.e. its electric field 
is perpendicular to the wavevector. Hence, it would seem inappropriate to use the present 
treatment. However, since the light wavevector is very small, one can think to rotate it so as 
to have it parallel to the electric field and apply the present formalism [65]. The validity of 
this approach has been rigorously demonstrated for cubic crystals in [13]. Furthermore, we 
observe that, in general, for anisotropic systems (3.3) and (3.4) depend on the direction of 
the vector q (i.e. on the polarization of the incoming radiation); thus, both microscopic and 
macroscopic responses are described by a dielectric tensor, instead of simple scalar functions. 

The same quantity €m is also related to electron energy-loss spectra (EELS) for vanishing 
momentum transfer, through the loss function — Im{l/€M}. For non- vanishing momentum 
transfer Q = q + G, the loss function is — Im{€ Q l G (q, co)}. In this case the longitudinal 
formulation of the dielectric response is obviously appropriate. From (3.3) it follows that 
the loss function can be related to the linear density-response function x • 

EEL(q + G,co) = -v G (q) Im {xcciq, co)} , (3.5) 

where the GG' -matrix x can be obtained by solving the Dyson-like screening equation (2.35). 
In (2.35) the full response function is expressed in terms of the independent-particle xks via 
a kernel composed of two terms: the bare Coulomb potential and the exchange-correlation 
contribution f xc . A similar expression can also be written in the case of the optical absorption, 
provided one builds a modified response function x ' 

Im{6 M } = - lim v G (q)lm{x GG }, (3.6) 

which satisfies the Dyson-like screening equation 

X = Xks + Xks(^ + Ac)x, (3.7) 
where the modified Coulomb interaction is defined as 
v G for 



v G 



(3.8) 

for G = 0. 



Following [14], the description of both absorption and EELS for q -> can be unified 
by introducing the generalized spectrum A(co) and a generalized function X GG (q,co). 
The function X stands for the modified response function x m the case of absorption and 



Time-dependent density-functional theory for extended systems 



371 



for the reducible response function x in the case of EEL: 




(3.9) 



In any formulation, the basic ingredients to obtain either the absorption or the EELS are the 
Kohn-Sham eigenf unctions and the eigenenergies that enter the expression for the independent- 
particle Kohn-Sham response function xks (2.30). These are usually obtained through a 
ground- state DFT calculation using an approximate exchange-correlation potential. In the 
total kernel of (2.35) and (3.7) v accounts for classical depolarization effects (also known as 
crystal local-field effects (LFE) in a solid). It reflects the microscopic-induced Hartree potential 
created by polarizable inhomogeneities in the system. The apparently subtle difference 
between the absorption and the EEL, i.e. the inclusion or exclusion of the long-range term 
l>o, is crucial for extended systems — for example, t>o is responsible for the plasmons, whereas 
its contribution in finite systems becomes vanishingly small [66]. The term f xc is a complex 
quantity that contains all non-trivial many-body effects. Its analytical expression is unknown. 

We have therefore two key approximations: (i) the ground- state exchange-correlation 
potential and (ii) the exchange-correlation kernel. Of course, these two quantities are in 
principle linked, due to the fact that the exchange-correlation kernel is the functional derivative 
of the time-dependent exchange-correlation potential. The relative importance of the two 
approximations depends, as we will see in the following, on the physical system under study. 
For example, when dealing with finite systems it is often essential to have good Kohn-Sham 
eigenstates — and therefore a good ground-state exchange-correlation potential — while the role 
of the exchange-correlation kernel is less relevant [22, 67] . The opposite is usually true for 
extended systems, where a good approximation to the exchange-correlation kernel turns out 
to be essential, particularly when it comes to describing optical absorption spectra [14]. 

When searching for approximations several approaches are possible. For the ground-state 
exchange-correlation potential, 40 years of development have led to a swarm of functionals 
(for more information, see one of the numerous review papers on the subject, e.g. [68]). For 
the exchange-correlation kernel, one can either look for a good approximation for the time- 
dependent exchange-correlation potential and then use the definition of f xc (2.34) or directly 
find an expression for the exchange-correlation kernel. Either choice has clear advantages 
and disadvantages. On the one hand, the exchange-correlation kernel is a simpler object, 
in the sense that it is a functional of the ground-state density, and is therefore amenable to 
more controllable approximations. On the other hand, if we are in possession of a good time- 
dependent exchange-correlation functional, we can tackle both linear and non-linear response 
properties. Instead, only linear response is accessible through the knowledge of the exchange- 
correlation kernel. Most often, the standard approximations used for v xc (r, t) are adiabatic 
and lead to very simple / xc s. 

3.1. Basic approximations: RPA and TDLDA 

The lowest level approximation to perform real calculations consists of setting to zero the 
terms in (2.35) and (3.7) coming from the microscopic components of the induced Hartree 
potential (v = 0) and the variations of the exchange-correlation potential (f xc = 0). By 
comparing (3.4) with (3.6), with vanishing v and / xc , it is possible to see that this is equivalent 
to neglecting all the off-diagonal components of the matrix e ~ G , . We will refer to this as the 
independent-particle approximation (IPA). The excitation energies in the IPA are simply given 
by the differences between the eigenenergies of the unoccupied and occupied Kohn-Sham 
states, which are used to build xks- Independently of the quality of the states entering in 
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Xks> this usually leads to absorption peaks that are systematically red-shifted in relation to the 
experimental spectra [14]. This is a consequence of the well-known gap problem of DFT: the 
gap between filled and empty states is substantially underestimated [69] . 

By neglecting only the exchange-correlation kernel we obtain the so-called random-phase 
approximation (RPA) [63,64]: 

A R C PA = 0. (3.10) 

In this case, the only part of the total kernel in (2.35) and (3.7) which is taken into account is 
the classical Coulomb term. This term describes the well-known Lindhard theory of screening 
with the addition of LFE [70] . Although very simple, the RPA yields results of reasonable 
accuracy for a wide range of systems, and it is still widely employed in actual calculations. 
We will see some examples in the following section. 

In the next step in the ladder of complexity come the already mentioned adiabatic 
approximations. At the level of the time-dependent exchange-correlation potential the 
adiabatic approximation implies 

„adiabatic W( ^ f) = j^r^)]^ ^ (3.H) 

where v xc [n] is some given ground-state exchange-correlation functional. Note that, regardless 
of the choice of v xc [n], the resulting kernel is instantaneous: / xc (f , r f , t, t f ) = 8(t — t f )f xc (r, r'), 
i.e. its Fourier transform is frequency independent (it can be of course non-local in space). If 
the LDA static potential is inserted in (3.11), one obtains the most common functional of 
TDDFT: the ALDA potential. Using definition (2.34), one can then derive the local and static 
ALDA (also called TDLDA) exchange-correlation kernel [17, 18]: 



j 2 HEG/ \ 

/ x T c DLDA (r, r', t, t') = 8(r - r')8(t - t f ) xc 1 j 



dn 2 



(3.12) 

=n G s(r) 



where e UEG is the energy per unit volume of the homogeneous electron gas and hqs * s the 
ground- state electron density. 

It is clear that the TDLDA retains all the problems already present in the LDA. The most 
important of these are perhaps, for neutral finite systems, the incorrect asymptotic behaviour 
of the LDA potential (instead of decaying as — 1/r, the LDA exchange-correlation potential 
goes to zero exponentially) and, for infinite systems, its local dependence on the density. These 
drawbacks cannot be corrected by using in (3.1 1) most of the GGAs [71] or the more modern 
meta-GGA functionals [72,73]. 

Nevertheless, for the calculation of optical response spectra in a large variety of finite 
systems, the TDLDA has proved to be able to reproduce low energy peaks with an accuracy 
of around 0. 1-0.4 eV [74]. For solids, the situation is a bit more complicated. EEL or x-ray 
scattering (IXS) spectra are often of good quality, particularly when the transferred momentum 
q is finite [25,26,75]. Instead, the description of Im{€ M (&0} for vanishing momentum transfer, 
which is the case in optical response, is perhaps the best-known failure of TDLDA. Note that 
in an extended system the TDLDA kernel for a vanishing q always yields a relatively small 
correction to the RPA results, because it is constant for q — ► and multiplied with xks that 
goes to as q 2 . It can hence only have an effect via the LFE [27]. We will discuss this issue 
in more detail in the following section. 

It is clear from the above that the behaviour of the different approximations depends 
strongly on the spectroscopy and on the dimensionality of the physical system. Therefore, 
it is interesting to present an overview of results for finite systems (molecules, clusters) to 
be compared with analogous results for infinite systems. In view of that, the following two 
subsections separately handle optical absorption in finite (section 3.2) and periodic systems 



Time-dependent density-functional theory for extended systems 



373 




Figure 2. Photo-absorption spectrum for a Nas cluster [14]. Dots: experimental data [82]; dashed 
line: IPA Kohn-Sham transitions; solid line: RPA; dashed-dotted line: TDLDA. 



(section 3.3). We then turn to EEL and IXS in section 3.4. Finally, we try to understand the 
failures of the TDLDA (section 3.5) and discuss possible routes to overcome them. 



3.2. Finite systems 

The first calculations of excitation energies within TDDFT were performed before the 
formal demonstration of the Runge-Gross theorem. In 1977 Ando determined intersubband 
transitions in semiconductor heterostructures [76]. Shortly after, Zangwill and Soven [15] 
applied TDLDA to the calculation of the photo-absorption cross section of rare gas atoms, 
obtaining very good agreement with experimental data. 

However, it was only in recent years that TDDFT became one of the most popular tools 
for the calculation of excitation properties. By now, the TDLDA kernel (3.12) has been 
successfully applied to atoms, organic and biological molecules, metallic and semiconducting 
clusters, fullerenes, etc [74,77-80]. Besides some more problematic cases (see section 3.5) the 
calculated excitation energies and absorption spectra are, in general, in excellent agreement 
with available experimental data. From the plethora of available applications of TDDFT to 
finite systems [81], we show here only two illustrative examples: (i) Nag, a prototype metallic 
cluster, and (ii) the poly cyclic aromatic hydrocarbons (PAHs). 

The absorption spectra of an Nag cluster [14] are shown in figure 2. The IPA calculation 
(dashed line) is compared with an RPA (solid line) and a TDLDA result (dashed-dotted line). 
The experimental spectrum [82] is also plotted (dots). As a general feature, present in both 
metallic and semiconducting clusters, we find that the RPA peaks are blueshifted with respect to 
the IPA peaks. Adding the TDLDA kernel brings a further correction, this time as a shift towards 
lower energies. The resulting transition energies then accurately reproduce the experiment. It is 
important to observe that, already in the RPA, absorption at low energies is correctly suppressed 
with respect to IP calculations. In fact, this suppression of the oscillator strength is essentially 
due to the induced classical depolarization potential. From (3.4) it can be observed that the 
LFE come from the off-diagonal terms of the matrix € GG > . In other words, they express the fact 
that the electronic response of an inhomogeneous structure is position-dependent (and not only 
distance-dependent). It is intuitive that such an effect has to be stronger the larger the inhomo- 
geneity of the system. Since a finite object represents a strong inhomogeneity in the otherwise 
empty space, it is not surprising that the LFE are particularly important for this kind of system. 
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Figure 3. Comparison between the calculated weighted sum of the spectra [83] (dotted-dashed 
line) for neutral PAHs and the experimental results of Joblin et al [84] (crosses) for a natural mixture 
of neutral PAHs with average Nq = 24. 

To illustrate the wide range of applicability of TDDFT, we now turn to the modelling of 
interstellar photophysics [83]. Due to their spectral properties, their high photo- stability and 
the fact that they are carbon-based, free gas-phase PAHs (in different charge and hydrogenation 
states) are commonly thought to be an important component of the interstellar medium. 
However, to compare with the measured spectrum of the interstellar medium, laboratory data 
of the individual PAHs are necessary. These are, however, often scarce and difficult to obtain. 
An alternative is numerical experiments based on TDDFT. 

In figure 3 the dotted-dashed line represents the calculation by Malloci et al [83] of 
the photo-absorption cross section up to the vacuum ultraviolet for a mixture of 20 different 
neutral PAHs, ranging in size from naphthalene (CioHg) to dicoronylene (C48H20). The overall 
spectrum for the PAH mixture is obtained as a weighted sum of the spectra for the single 
C w H m molecules. (The statistical weights are assumed to be inversely proportional to the total 
number of carbon atoms Nq of each molecule. The average Nq of the theoretical sample is 
Nq = 23.55.) The experimental spectrum measured by Joblin et al [84] for a natural mixture 
of neutral PAHs with average Nq = 24 is also plotted in figure 3 (crosses). Two distinct 
features can be observed in the spectrum: (i) a collective broad absorption peak, resulting from 
the sum of the tc -> n* transitions at ~6 eV and characterized by distinct structures due to the 
coincidence of relatively strong transitions in different molecules, and (ii) a smooth far-UV 
rise. The agreement between theory and experiment is very good and validates the use of 
TDLDA calculations as a substitute for laboratory data when the latter are lacking. 

We have just witnessed snapshots of the quality of TDLDA calculations for finite systems. 
However, it should not be forgotten that in some cases TDLDA is not adequate to describe 
excitations of finite systems: a typical example is the failure in reproducing Rydberg series [22] . 
Moreover, when the molecules become more extended, this quality in general degrades. An 
example is the calculation of optical properties of long conjugated polymers [24, 85]. The 
problem is related to a non-local dependence of the exchange-correlation potential: in a system 
with an applied electric field, the exact exchange-correlation potential develops a linear part 
that counteracts the applied field [24, 86]. This term is completely absent in both the LDA and 
the GGA. It is present in more non-local functionals such as the EXX (see section 7). 
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Figure 4. Imaginary part of the macroscopic dielectric function for LiF [87] . Dots: experiment [88] ; 
dotted line: BSE calculation; dashed line: IPA calculation; solid line: RPA calculation; dotted- 
dashed line: TDLDA calculation. 

In the following we will investigate this problem more closely. 

3.3. Optical absorption in extended systems 

The simplest approach to the optical properties of semiconductors or wide-gap insulators 
within TDDFT is the TDLDA. In view of the excellent quality of the results obtained for the 
photo-absorption of clusters, one could perhaps expect that the same would occur for extended 
systems. This is unfortunately not the case. As we can see in figure 4 for the optical absorption 
of a LiF crystal, the TDLDA (dashed-dotted line) induces only some minor modifications with 
respect to the RPA (solid line), and both are very far from the experimental curve (dots). The 
largest disagreement concerns the absence of the strong excitonic peak at about 12.5 eV. Good 
agreement can be found using the many -body Bethe-Salpeter approach (BSE, dotted line) 
at the price of significantly larger computational effort [89-92]: in that framework, electron 
addition and removal energies as well as the electron-hole interaction are explicitly calculated 
within many-body Green's function theory (see section 7). 

This situation is quite general and is found in a wide range of semiconductors (Si — see 
also the inset of figure 5 — Ge, GaAs, etc) and wide-band gap semiconductors or insulators 
(diamond, MgO, S1O2, etc). It is typical for absorption, as opposed to loss spectroscopies, 
even when both techniques are employed to study the same system. A detailed analysis of the 
problem will be the subject of section 6.2. 

3.4. EELS and IXSS in extended systems 

In figure 5 we can observe both the absorption and EELS at vanishing q (within RPA) for bulk 
silicon [66]. To interpret this picture it is useful to use the generalized spectrum A(co) of (3.9). 
The modified RPA polarization function X (3.9) can be further generalized as 

X(co) = (1 - XksK^o - Xks^) _1 Xks- (3.13) 

If y = 1, A(co) = EELS, and if y = 0, A(co) = Abs. Moreover, it is possible to follow the 
evolution of the spectrum when y varies continuously from 1 to 0. Figure 5 shows how the 
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Figure 5. Continuous connection between EELS and absorption spectrum of bulk silicon, via 
vo [66]. Experiments from [93] (absorption) and [94] (EELS). 

EELS turns continuously into the absorption when vq is switched off. This exemplifies the 
action of the long-range component vq, which is responsible for the huge difference between 
an EELS and an absorption spectrum. 

Let us go back to the EEL spectrum of bulk silicon: a comparison between the experimental 
and the RPA spectrum is shown in the second inset of figure 5. Olevano and Reining [14, 95] 
showed that the TDLDA has better agreement with experiment than RPA, even though the 
difference is small (see figure 16 of [14]). Some improvement can be found when the BSE 
approach is used (see again figure 16 of [14]). But since the full BSE calculation of a valence 
plasmon is still a computationally involved task, the use of TDLDA (or even RPA) is often 
well justified. 

When the electron density does not present particular inhomogeneities, it can be enough 
to include only vq in the kernel (2.35) to obtain an accurate calculation of the loss spectra 
of extended systems. In the case of layered or low-dimensional structures [96], or in the 
presence of localized states [97, 98], the contribution of v also becomes essential and only 
a RPA, or often better TDLDA, calculation can yield good agreement with the experimental 
data. The similarity between RPA and TDLDA is a quite general feature for the loss function 
at small transferred momentum q. It holds, for example, for the loss function of graphite [96] 
and for the integrated loss function of TiC>2 [97]. Good agreement with experimental spectra 
was obtained within RPA also for the EELS of diamond [27] and ZrC>2 [98], always at low 
momentum transfer. As a general rule, when q gets larger, the contributions of LFE and of 
the exchange-correlation kernel within TDLDA become more important. Also in this case 
RPA and TDLDA allow good agreement with experiment, as is shown by recent calculations 
of IXSS at the RPA and TDLDA level for Al [25], rutile Ti0 2 [75] and various 3d transition 
metals [26]. In some cases, TDLDA can give a sizable improvement with respect to RPA [28]. 
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Figure 6. Integrated loss function of T1O2 for (q ~ 0.4 A -1 ) [97]. Solid line: experiment. Dashed 
line: RPA calculation. The inset compares the IPA calculations (dashed line) with experiment (solid 
line). 



The case of T1O2, illustrated in figure 6, is particularly interesting: at small momentum 
transfer the IPA spectrum yields results up to about 25 eV very similar to those obtained 
within the RPA or the TDLDA. At higher energies the IPA picture breaks down as it cannot 
describe correctly the structure originated by transitions from the localized semi-core states. 
The agreement with the experimental data is obtained only within RPA and TDLDA (the 
TDLDA spectrum is very similar to the RPA one and therefore it is not shown in figure 6). 

These statements often hold when the momentum transfer is small. For large momentum 
transfer, the LFE become more important and, with the increasing contribution of G ^ 
terms also the influence of the TDLDA kernel increases. In fact, by comparing IPA, RPA and 
TDLDA calculations of IXSS for T1O2, Gurtubay et al [75] proved that at large momentum 
transfer the calculations agree with the experiment only when LFE are included, even at low 
energies (see figure 7). 

In conclusion, we can state that the TDLDA is often very reliable for EELS and IXSS 
(both for small and large momentum transfer) and for photo-absorption in finite systems. LFE 
often give a sizable contribution to this success. One of the main remaining problems is the 
optical absorption in extended systems. 



3.5. What is missing in RPA and TDLDA? 

We can summarize now the situation as it was explored so far: 

• For excitation properties of finite systems, in general, RPA and TDLDA work quite well. 
There are of course many exceptions, most of which are related to the incorrect tail of 
the LDA (or GGA) exchange-correlation potential at large r. Some problems related 
to this deficiency are the already mentioned impossibility to reproduce Rydberg series, 
the overestimation of polarizabilities in long chain molecules, the large underestimation 
of ionization energies or the wrong description of any situation where the electrons are 
pushed to regions far away from the nuclei (e.g. by a strong laser). These issues can be 
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Figure 7. Dynamical structure factor of T1O2 at \q\ =1.32 A -1 along the [00 1] direction [75]. 
Solid (dotted) line: TDLDA spectrum with (without) LFE. Dashed line: RPA spectrum with LFE. 
Circles: IXS measurements [75]. 



solved by the use of functionals with the correct asymptotic behaviour, such as the EXX 
or the adiabatic LB94 [99]. 

• For extended systems, EELS and IXSS at small and large momentum transfer are often well 
reproduced within TDLDA. Instead, TDLDA fails in the calculation of optical (q = 0) 
spectra of non-metallic solids [23] . To explain this failure, the wrong asymptotic behaviour 
of the exchange-correlation potential is less relevant, while the wrong asymptotic limit 
of the exchange-correlation kernel is crucial. For infinite systems, the q = component 
of Xks vanishes as q 2 . It is then clear from the response equation (3.7) that if f xc has to 
correct the non-interacting response for q -> it will have to contain a term that behaves 
asymptotically as l/q 2 when q — >► 0. This term will be particularly important in absorption 
calculations, where the Coulomb part of the kernel does not contain vo ~ 1/g 2 . This 
crucial term cannot be found in the local-density or gradient-corrected approximations. 

Finally, we recall that the TDLDA exchange-correlation potential is local in time. Few 
attempts to derive functionals which are non-local in time, i.e. which include memory effects, 
have been made so far. By analogy with hydrodynamics, Dobson et al assumed that in the 
electron liquid memory resides not with each fixed point r, but rather within each separate 
'fluid element' [100]. Thus, the element which arrives at location r at time t 'remembers' 
what happened to it at earlier times when it was at locations different from its present location 
r. Using this concept, Dobson et al proposed a functional that satisfies Galilean invariance 
and Ehrenfest's theorem. Unfortunately, no applications of this functional exist to date. 
This approach was further extended by Tokatly within time-dependent current DFT [101]. 
Furthermore, the frequency dependence of the exchange-correlation kernel has been proved 
to be essential to describe charge transfer between open-shell species [102] and double 
excitations [103-107]. An example of a frequency-dependent model exchange-correlation 
kernel will be presented in section 8.1. 

Several attempts have been made to correct the shortcomings of RPA and TDLDA. In the 
next section we will start by considering the case of metallic systems and discuss an explicit 
density functional beyond the ALDA. 
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4. Explicit density functionals 

Explicit density functionals are expressions for the exchange-correlation kernel which are 
defined directly in terms of the electron density. Most of the commonly used approximations, 
such as the ALDA, belong to this class of functionals. The majority is derived from numerical 
results for the homogeneous electron gas which are cast into a parametrized form, thus allowing 
a transfer to other systems. The prevalence of explicit density functionals is largely due to the 
fact that their evaluation is typically cheap and adds little computational overhead to practical 
calculations. However, the physical content of a given parametrization, particularly when 
transferred to a very different material, and its accuracy for the study of excited states are 
not always clear. Some problems of the RPA and the ALDA were already discussed in the 
previous section. In the following we consider a wider range of explicit density functionals that 
go beyond these basic approximations and analyse their performance for extended systems. 



4.1. Dynamic exchange-correlation effects in the electron gas 

In the following we examine the influence of the exchange-correlation kernel on the excitation 
spectrum of the homogeneous electron gas, based on [108]. As the effective ground-state 
potential is a trivial constant for this model system, the choice of the kernel is the only source 
of errors, whose impact can thus be clearly identified. In addition, the relative simplicity of the 
homogeneous electron gas makes it possible to explore accurate numerical constructions that 
go significantly beyond basic approximations such as the RPA and ALDA. Note however that 
their applicability to inhomogeneous semiconductors with a finite band gap, with which most 
of this review is concerned, is an entirely different question, because the electronic properties 
of metallic and non-metallic materials deviate fundamentally. Specifically, we here consider 
the following schemes. 

(a) In the RPA (3.10) all dynamic exchange and correlation effects are ignored by setting the 
kernel to zero. 

(b) The ALDA replaces the wavevector- and frequency-dependent kernel of the electron gas 
by its long wavelength and static limit (for inhomogeneous systems, this value is then 
used at each point in space according to the local density, see (3.12)): 

/ x f DA = lim f» EG ( ? ,« = 0). (4.1) 

Note that for the homogeneous gas the exact kernel in reciprocal space is of course not 
a matrix but has only a scalar dependence on the absolute value of q. Also note that, 
contrary to the case of semiconductors and insulators mentioned in earlier sections, no 
l/q 2 divergence appears. 

(c) In their original application of TDDFT to excited states Petersilka, Gossmann and Gross 
(PGG) [20] derived an exchange-only kernel within the approximation of Krieger et al 
[109]. Designed for small atoms, the PGG formula is, in fact, exact for two-electron 
exchange, but deviations are expected for extended systems. In particular, it does not 
contain the frequency dependence of the exact exchange kernel [110]. The formula for 
the non-local PGG kernel f^ G (q) in reciprocal space for the homogeneous electron gas 
is given in [111]. 

(d) Burke, Petersilka and Gross (BPG) [112] proposed a hybrid formula that further improves 
the excitation spectra of atoms by incorporating correlation as well as a self-interaction 
correction. It combines expressions for symmetric and antisymmetric spin orientations 
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from different approximations in a spin density-functional formalism. For the unpolarized 
homogeneous electron gas this kernel reduces to 

A B c PG (<7) = + f£ftl (4-2) 

(e) A good parametrization of the static exchange-correlation kernel for the homogeneous 
electron gas was given by Corradini, Del Sole, Onida and Palummo [113] (CDOP), 
who used the Monte Carlo results of Moroni et al [114] for the static local-field factor 
G(q) = — f^ c (q)/v{q). Unlike the original data, this parametrization is not restricted to 
metallic densities, because it incorporates the known asymptotic limits for high and low 
densities. By construction, the CDOP kernel becomes identical to the ALDA in the long 
wavelength limit. 

(f) Finally, we consider a parametrization of the dynamic local-field factor of the 
homogeneous electron gas proposed by Richardson and Ashcroft [115] (RA), including 
the corrections given in [1 1 1], which stems from the summation of self-energy, exchange 
and fluctuation terms in the diagrammatic expansion of the polarization function. It 
satisfies many important sum rules and reproduces the exact asymptotic expressions for 
small and large wavevectors. At intermediate wavevectors and frequencies it provides a 
realistic description of the position and magnitude of extrema, which are related to the 
pair distribution function evaluated at zero separation. Because of this careful derivation 
one can expect the RA expression to be very close to the exact dynamic exchange- 
correlation kernel of the homogeneous electron gas and to give an accurate account of 
the plasmon dispersion. In the absence of experimental data, we therefore use the RA 
results as a reference in order to assess the performance of simpler approximations. The 
parametrization was originally given on the imaginary frequency axis; here we use its 
continuation to the full complex plane. 

Among the many other expressions of practical relevance we mention the early 
parametrizations of the static local-field factor by Hubbard [116], Vashishta and Singwi 
[117, 118] and Utsumi and Ichimaru [119], which are now superseded by more accurate 
Monte Carlo results and an attempt by Gross and Kohn [17] to retain the frequency dependence 
within the local-density approximation at long wavelengths. Unfortunately, the latter violates 
a number of exact conditions, particularly the harmonic-potential theorem [120] . The dynamic 
and non-local exact exchange kernel, an implicit density functional constructed from the 
Kohn-Sham wavef unctions, as well as other orbital-dependent functionals, are discussed later 
in section 5. 

The poles of the linear density -response function of the homogeneous electron gas 

X(q,co) = (4.3) 

1 - Xks(<7, co)[v(q) + f xc (q, co)] 

stem from two different sources. The singularities of the Kohn-Sham density-response 
function in the numerator correspond to independent electron-hole pair excitations; they form 
a continuum that is bounded by the lines \q 2 — qk? ^ co ^ \q 2 + qk?. In addition, the zeros of 
the denominator give rise to a distinct plasmon branch co q , which describes resonant collective 
charge oscillations of the electron system. We focus on the latter, since the plasmon dispersion 
gives a direct measure for the quality of the kernel. Before scrutinizing the numerical results 
we first discuss what can be deduced from an analytic expansion of the plasmon dispersion up 
to second order in q [121]: 



: CO p ] 
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(4.4) 
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Figure 8. Plasmon dispersion for the homogeneous electron gas at r s = 4 calculated with different 
approximations for the exchange-correlation kernel (see text). The grey-shaded area marks the 
electron-hole pair continuum. The finite imaginary part of the plasmon energy in this region is 
also shown (see also [108]). 



where &tf = 2(3ft/jr) 1/6 is the Thomas-Fermi wavevector. As none of the above f xc diverges, 
all curves approach the classical plasma frequency co p \ = (A7tn) l/2 in the long wavelength 
limit. The kernel only introduces corrections in second order, where the element f xc (0, co v \) 
appears. The ALDA contains by construction (4.1) the correct long wavelength limit, but 
its neglect of the frequency dependence introduces an error in the parabolic term. The 
CDOP formula produces the same second-order term as the ALDA since the two kernels, 
which are both static, coincide for lim^o- Also the PGG and BPG functionals are static 
approximations; moreover, they do not approach the correct long wavelength limit of the 
homogeneous electron gas and therefore generate a different parabolic coefficient. The RA 
kernel, which incorporates the full frequency dependence, is the only parametrization that is 
formally exact beyond the trivial zeroth order. 

The calculated plasmon dispersions for r s = 4, obtained from a numerical search for the 
zeros of the denominator of (4.3) in the complex frequency plane, are shown in figure 8. The 
results are representative of the whole range of metallic densities [108]. As predicted, all 
curves start at the classical plasma frequency. For small wavevectors only a minor spread of 
the results is observed, because the factor 9/(10^ F ) in (4.4) outweighs the contribution of 
the kernel. However, a slight downward shift compared with the RPA is clearly visible for 
all non-trivial approximations, because dynamic exchange and correlation effects combine to 
lower the energy of the electron system. The ALDA and the CDOP formulae produce curves 
that are initially very close to the reference RA result, indicating that the neglected frequency 
dependence is of little consequence as long as the correct long wavelength limit is reproduced. 
This point is emphasized further by the relatively large deviation for the static PGG kernel, 
which stems precisely from its incorrect behaviour at q —> 0. The BPG curve, as expected, 
lies between the ALDA and PGG results. 

At larger wavevectors, where the parabolic expansion (4.4) is no longer valid, the 
differences between the considered approximations become more pronounced. The dispersion 
resulting from the static ALDA kernel remains close to the RPA and yields too high energies, 
while the CDOP result begins to deviate slightly from the RA curve after the onset of damping 
in the electron-hole pair continuum. This discrepancy must be attributed to the static nature of 
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the CDOP kernel. Furthermore, it can be seen that the strong downward shift of the exchange- 
only PGG formula leads to an even larger error in absolute terms than the underestimation 
of dynamic exchange and correlation effects in both the ALDA and the RPA. The hybrid 
BPG formula, which combines the PGG and ALDA parametrizations, profits from a partial 
cancellation of errors but improves only marginally upon PGG. 

Due to the possible decay into electron-hole pairs, the plasmon energy contains a non- 
zero imaginary part, also displayed in figure 8, whose inverse is proportional to the excitation 
lifetime. As a general rule, all kernels yield the same quality of approximation for the imaginary 
part as for the real part of the plasmon energy. At small wavevectors all static kernels predict 
a vanishing imaginary part, which corresponds to an unphysical infinite lifetime. This artefact 
results from modelling f xc (q , co) as a purely real quantity by evaluating it at co = 0. In contrast, 
the exact kernel has a finite imaginary part at non-zero frequencies, which for small wavevectors 
is related to the multi-pair component of the linear density-response function [122]. Such 
multi-pair decay channels are ignored in the RPA and related schemes, which is ultimately 
the reason for their qualitatively wrong behaviour. Mermin's modification of the Lindhard 
dielectric function avoids the problem of infinite lifetimes [123], but the correction based 
on relaxation times is introduced in a phenomenological manner that makes it unsuitable 
for ab initio calculations. Among the expressions considered here, only the dynamic RA 
parametrization correctly predicts a finite plasmon lifetime over the entire frequency range, 
although outside the electron-hole pair continuum the imaginary part of the plasmon energy 
is several orders of magnitude smaller than the real part and not discernible in the figure. 

The good agreement between the static CDOP parametrization on the one hand and 
the dynamic RA result on the other hand over a large wavevector and density interval 
indicates that the frequency dependence of the kernel plays a weak role for the plasmon 
dispersion in the homogeneous electron gas. In contrast, the significant discrepancy between 
static approximations such as the ALDA that contain the correct long wavelength limit and 
others, such as PGG, which do not, suggests that a correct parametrization of the wavevector 
dependence is crucial. Similar conclusions concerning the relative importance of the frequency 
and wavevector dependence were also reported for the correlation energy of the homogeneous 
electron gas [111]. 

4.2. Extension to inhomogeneous systems 

As the homogeneous electron gas allows a highly accurate treatment of the kernel that includes 
the full wavevector and frequency dependence, it is tempting to transfer these results to real 
inhomogeneous systems. Such an approach is mathematically justified if the variation of the 
density ftGs(f) = n + An(r) with f An(r) d 3 r = from a homogeneous charge distribution n 
is small, i.e. \An(r)/n\ <^ 1. In this case the kernel can safely be approximated as [3] 

/ xc (r, r\ co) * f* EG (n, \r - r'|, co). (4.5) 

It should be noted that rapid density oscillations leading to large gradients are not explicitly 
excluded, as long as the magnitude of the oscillations themselves remains sufficiently small. 

Unfortunately, the approximation (4.5) is of rather little practical use, because the condition 
of an almost constant electron density is almost never fulfilled for real materials. However, 
for co = a very similar expression with a much wider range of validity can also be obtained 
under less restrictive assumptions. The derivation is based on the observation that the static 
exchange-correlation kernel, although non-local, is actually short-ranged and decays rapidly 
if the separation \r — r f \ significantly exceeds the inverse Fermi wavevector. Furthermore, 
it depends only on the electronic structure in the vicinity of the points r and r' . This is 
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a manifestation of the 'nearsightedness' principle [124] and is caused by the destructive 
interference of wavefunctions in quantum-mechanical many-particle systems over a typical 
de Broglie wavelength. The nearsightedness principle is strictly valid only in equilibrium and 
hence does not apply to the dynamic exchange-correlation kernel at non-zero frequencies. 
Therefore, it fails to cover spectroscopies such as EELS or optical absorption that include a 
finite energy transfer, but the static kernel may be useful in other situations, such as total-energy 
calculations [3, 125]. 

In combination with the short range of the static exchange-correlation kernel, the 
nearsightedness principle implies that the approximation [3, 125] 

/ xc (r, r', 0) « /™ G (n(r, r'), \r - r'|, 0) (4.6) 

is valid if the density variation is small on a length scale given by a few inverse Fermi 
wavevectors. The definition of the average density n(r,r f ) is arbitrary from a mathematical 
point of view; as long as it reduces to h = hqs (?) = ^gs(^) in the limit of a homogeneous 
density distribution, the expression (4.6) is exact to the zeroth order. It is reasonable to 
expect, however, that a suitable choice might improve the performance of the functional 
for realistic inhomogeneities. Indeed, it has been demonstrated that some possible and, at 
first sight, plausible definitions such as the mid-point density n(r, r f ) = nQs(\(r + r f )) can 
even lead to divergences in total-energy calculations for strongly inhomogeneous systems, 
such as atoms and surfaces [126]. The reason for the observed unphysical behaviour is 
that h (r, r r ) approaches zero in this case if one spatial argument moves far into the vacuum 
region while the other remains inside the system. As a consequence, the corresponding 
Fermi wavevector also approaches zero, and the kernel becomes long-ranged, so that spatial 
integrations pick up many incorrect contributions. The problem is avoided by the choice 
n(r,r f ) = |[^gsW + ftGs(f')L which guarantees a finite density and Fermi wavevector in 
the same situation and is thus recommended. It can indeed be argued to constitute the most 
logical choice, because it represents a smooth function by its value at some average point of 
its arguments (WeierstraB) [125]. 

Finally, it should be borne in mind that the approximation (4.6) is still derived by an 
expansion starting from the homogeneous electron gas. Therefore, its validity is restricted to 
systems where perturbation theory is applicable, i.e. where the inhomogeneity of the density 
does not alter the physical properties of the system qualitatively. This includes, for example, 
many bulk metals. On the other hand, it fails to improve the optical absorption spectra of 
semiconductors, such as silicon [127], because the Fourier transform of (4.6) tends to the 
finite value of the ALDA kernel in the long wavelength limit. Strictly speaking, of course, the 
restriction of the nearsightedness principle to static phenomena already precludes the use of 
(4.6) from the outset in this case. 

5. Orbital-dependent functionals 

In contrast to explicit density functionals, orbital-dependent functionals are constructed from 
the Kohn-Sham wavefunctions (of course, they still depend implicitly on the electron density 
through the self-consistency condition). Although the computational cost is typically much 
higher than the straightforward evaluation of a parametrization in terms of the density, this 
approach has the advantage that it offers a systematic route to successively more accurate 
approximations. Orbital-dependent functionals can be obtained in various ways. One can get 
a potential and a kernel from a functional derivative of a suitable approximation of the action, 
for example, by expanding the latter in powers of the Coulomb interaction; this is outlined 
in section 5.2. Alternatively, one can exploit links between many-body perturbation theory 
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(MBPT) and the density-functional formulation. In the next subsection we briefly show how 
known results can be obtained from a linearized version of this link; section 6 is dedicated to a 
more detailed discussion of this connection, including results that do not involve a linearization. 



5.7. Potentials and kernels from a linearized Sham-Schliiter equation 

The density-functional and the MBPT framework are linked by the requirement that the former 
has to yield the correct density given by the one-particle Green's function G of the latter, via 
n(r, t) = — iG(r, r f ,t,t + ), for the static as well as for the time-dependent case. (Note that 
in the latter case it is recommended to use Keldysh Green's functions in order to obtain the 
physical densities via the simple relation above [128].) From the Dyson equation one therefore 
obtains 

= Zd2d3 G KS (1, 2) (S(2, 3) - v xc (2)8(2, 3)) G(3, 1 + ), (5.1) 



where the number '1' represents position, spin and time (r\, &\,t\) and £ is the self-energy 
[11, 129, 130]. This is the so-called Sham-Schluter equation [131, 128]. If G is replaced by 
Gks everywhere (including in the construction of E, as symbolized by S KS ) the solution of 
this linearized equation yields for the potential 



^xc(l) 



-i jd2 d3 j d4 x K s(^ 2)G KS (2, 3)E ks (3, 4)G KS (4, 2 + ). (5.2) 

The kernel / xc , the functional derivative with respect to the density of v xc , has a contribution 
/ (2) that stems from the derivative of E KS and other terms coming from the derivative 
of the Green's functions and the inverse response function [132]. Both are explicitly orbital- 
dependent. 

As an example, one can use the simplest case where E KS is approximated by the (KS) 
Fock operator E^ s (l, 2) = iGics(l, 2)i>(2, 1). In that case, (5.2) yields the so-called exact 
exchange (EXX) OEP potential [110, 133]: 

ifcxx(l) = "i fd2 d3 j d4 Xks(1> 2)G ks (2, 3)E x ks (3, 4)G ks (4, 2 + ), (5.3) 

where all KS quantities are calculated self-consistently using the EXX potential. The 
contributions to the kernel become 

/( 2),exx (1? 2) = j d3 d4 d5 d6xKs i ( i 9 3)G KS (3, 4)G ks (5, 3) 

xv(4, 5)G KS (4, 6)G ks (6, 5) Xks 1 (6, 2) (5.4) 

(which is nothing else but the electron-hole attraction term of time-dependent-EXX, TD-EXX 
[1 10, 133, 132]). The rest of the terms — which has the difficult task to open the band gap with 
respect to the KS one in TD-EXX — reads 

/ (D,exx (1) 2) = J d3 d4 d5 d6 x -i (lj 3)Gks(3; 6)G ks (6, 4) 

x[Sf (4, 5) - 5(4, 5)i; K s(4)]Gks(5, 3)/ks(6, 2) 



■jd3 d4 d5 d6 Xks(1, 3)G ks (6, 3)G ks (3, 4) 



x[Ef (4, 5) - 5(4, 5>ks(4)]G K s(5, 6) X ^(6, 2). (5.5) 
TD-EXX will be discussed more in detail in the next subsection. 
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5.2. The time -dependent optimized-potential method from an action formalism 

The TDOEP potential can also be obtained from the action formalism. In this case, the time- 
dependent optimized-potential method [49] treats the Coulomb interaction as a perturbation 
that is switched on adiabatically in the interval (— oo, to), while a compensating local potential 
ensures that the density remains constant and equal to the static ground-state density ncs (/*) 
throughout the entire switching-on process. In this way it provides an adiabatic connection 
between the stationary ground state of the non-interacting Kohn-Sham system at t -> — oo and 
the wavefunction ^ of the true interacting electrons at t = to that enters in the definition of the 
action (2.21). In order to incorporate the switching-on process in the theoretical description 
of the evolution of the system, the beginning and end of the pseudotime contour C must be 
extended to — oo. The combination of the adiabatic connection with the time-contour method 
makes it possible to apply standard perturbation techniques and expand A xc [n] in terms of the 
Kohn-Sham orbitals and the Coulomb interaction [128]. 

For an orbital-dependent action, two equivalent expressions for 8A xc [n]/8v^s(r, x) can 
be derived by applying the chain rule with different intermediate quantities: 



jf drV(r') j 



, 3 f 8A xc [n] 8n(r f ,x f ) 

a r 

8n(r f , x')8v KS (r, x) 

= j drV(r') fdVf;f^^ + 8Axcl(fj] ^ {r '' X ' y 




x') 8v KS (r, x) 8cp*(r>, x') 8v KS (r, r) ) 

(5.6) 

The first functional derivative on the left-hand side gives rise to the exchange-correlation 
potential (2.25), which can be determined from this identity because all other terms are known: 
the linear Kohn-Sham density-response function 8n(r f , x f ) /8vjts(r, x) is obtained with the help 
of time-dependent perturbation theory by applying the same techniques as in section 2 in the 
pseudotime domain, while the derivatives of A xc [ft] on the right-hand side of (5.6) can be 
calculated analytically for a given orbital-dependent functional. The variations of the orbitals 

S f i(r / T l = ~ i( Pj( r > *) y>j?(r, *)n(r\ r')0(r' - x) (5.7) 
8v KS (r, x) 

and the corresponding formulae for the conjugate orbitals 8cp*(r f , x f )/8vKs(r, x) again follow 
from time-dependent perturbation theory. Finally, an expression for the exchange-correlation 
kernel is obtained by manipulating in an analogous way the second functional derivative 

8 2 A xc [n] f f 3 3 8 2 A xc [n] 8n(r\ rQ 

= / art {x )ax\t (rj / d r a n 

8v KS (r, T)8v KS (r 2 , r 2 ) J c J 8n(r' , x')8n(n, n) 8v KS (r, x) 

8n(r u xi) f f 3 , 8A xc [n] 8 2 n(r f ,x') 



■ j &x't'(x') j d 3 r 



8v KS (r 2 ,x 2 ) Jc J 8n(r f , x') 8v KS (r, x)8v KS (r 2 , x 2 ) 

(5.8) 

here the term on the right side containing the double derivative of A xc [/t] is nothing else but 
Xks/xcXks> which allows one to solve for f xc after the evaluation of all other terms. 

5.3. Exact exchange 

The leading term in the expansion of A xc [ft] in powers of the Coulomb interaction is the 
exchange part, which is of first order and given by 

AAn] = ~2j c dxt ^ L fj 2^ fk ) d r J d r \^7\ ' (5 - 9) 

j = l k= 1 
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while the correlation part A c [n] includes all higher-order contributions. Inserting this 
expression into (5.6) and evaluating all quantities with the static density ftGs(f) on the physical 
time axis yields the 'exact exchange potential' (5.3) [134, 135]. For the static ground-state 
potential this yields a relatively simple expression: 



/ 



dV Vexk(/)xks(/, r; co = 0) 

-/*V kxw^S/W,)?^^. (5..0, 

7 = 1 k=l S -> Sk 

with the non-local exchange self-energy 

\r — r'\ 

As pointed out above the latter takes the same form as the Hartree-Fock exchange operator 
but is constructed here from the Kohn-Sham orbitals. Solving (5.10) for the exact exchange 
potential requires an inversion of the static linear density-response function xks ( r \ r \ 0), which 
is in practice achieved by a matrix inversion after projection on a suitable basis set, such as 
plane waves for extended systems [136, 137]. The matrix elements with reciprocal lattice 
vectors G = or G' = and q = must be omitted in this case [138]. The restriction to 
this submatrix is necessary to guarantee a one-to-one correspondence between variations of 
the density and the effective potential and hence to obtain an invertible operator; it excludes 
constant potential shifts that leave the density unchanged as well as density variations that 
violate particle-number conservation. 

Compared with the evaluation of explicit density functionals, the construction of uexxW 
is considerably more expensive, because it requires not only the occupied conduction states 
but also a summation over the unoccupied part of the spectrum. In addition, the linear density- 
response function must be inverted in each cycle of the self-consistency loop. Practical 
calculations were therefore initially restricted to atoms and small molecules, but the method 
is now also routinely applied to bulk semiconductors, insulators and metals [138-144]. 

From a theoretical point of view, the exact exchange potential has definite advantages. In 
particular, it is free of self-interaction and exhibits the correct asymptotic behaviour outside 
finite systems [135], where it decays like — 1/r , while the LDA and GGA fall off exponentially. 
The exact exchange potential further features a discontinuity with respect to a change in the 
number of electrons [143]. Such a discontinuity is also contained in the exact functional; it is 
the difference between the Kohn-Sham eigenvalue gap and the true quasiparticle band gap in 
semiconductors [145, 131]. 

The 'exact exchange' kernel is similarly obtained, following the outline of (5.8), from the 
relation 

f dV f dV'xKsOr, r"\ co)f x (r", #"'; co) X Ks(r"\ r'; co) = Pi(r, r'; co). (5.12) 

The right-hand side equals the first-order contribution to the irreducible polarizability in an 
expansion in powers of the Coulomb interaction. It consists of five distinct terms, which can 
be derived in matrix notation [110] but are more easily summarized in terms of Feynman 
diagrams [146]. The diagrammatic form of P\(r, r'\ co) is shown in figure 9 together with the 
representation of equation (5.10). The first two terms are the self-energy insertions with the 
non-local exchange operator for independent electrons (contributions due to E x in (5.5)), while 
the third arises from the attractive electron-hole interaction (5.4). The last two terms contain 
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Figure 9. Diagrammatic form of (a) the equation for the exact exchange potential and (b) the first- 
order irreducible polarizability. A solid line represents the Green function of the non-interacting 
Kohn-Sham electrons, a broken line represents the Coulomb interaction and the local exchange 
potential is indicated by a cross. 




Figure 10. Plasmon dispersion for the homogeneous electron gas at r s =4 with the exact exchange 
kernel (EXX) compared with the ALDA and the reference values obtained with the parametrization 
by Richardson and Ashcroft (RA) [115]. 



the exact local exchange potential (contributions due to v xc in (5.5)) and must be subtracted in 
order to avoid double counting. 

In contrast to the ALDA, the exact exchange kernel is both non-local and frequency- 
dependent. Furthermore, in the case of semiconductors its Fourier transform diverges in the 
limit of small wavevectors, as required for the exact functional [110]. 

As a first practical test of its performance for electronic excitations in extended systems 
we display the plasmon dispersion for the homogeneous electron gas at r s =4 calculated 
with the exact exchange kernel in figure 10. The ALDA and the results obtained with the 
parametrization by Richardson and Ashcroft [115] are shown for comparison and are the same 
as in figure 8. The exact exchange kernel is evidently in very good agreement with the RA 
reference values in the entire region outside the electron-hole pair continuum and constitutes 
a definite improvement over the ALDA. Based on our earlier analysis, we attribute this to a 
close match with the wavevector dependence of the complete functional and conclude that the 
exact exchange kernel is a suitable starting point for the quantitative investigation of collective 
excitations in free-electron metals [147]. 

Early results from the time-dependent exact exchange method for semiconductors 
indicated a good performance: Kim and Gorling [133] calculated the optical absorption 
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spectrum of silicon and found good agreement with experimental data. However, it was later 
observed that these calculations contained a cutoff of the divergent Coulomb potential at small 
wavevectors, which actually had a large effect on the results and was in fact responsible for the 
good quantitative performance [132] . As the singularity reflects the long range of the Coulomb 
interaction in real space, this cutoff means that the kernel was effectively evaluated with a short- 
range, screened interaction. If the singularity is properly taken into account, TD-EXX should 
indeed be similar to time-dependent Hartree-Fock. The results then change drastically, which 
can be understood by considering the effect on the polarizability P\ (r, r'\ co): the Kohn-Sham 
eigenvalues obtained from the exact exchange potential are closer to the experimental band 
structures than to the LDA ones [138], but the self-energy insertion in the first two terms in 
figure 9 increases the eigenvalue gap close to the Hartree-Fock result, leading to a far too 
high absorption threshold. In addition, the unscreened electron-hole interaction in the third 
term gives rise to a strongly overbound exciton with an incorrect line shape [132]. If the 
bare Coulomb potential is replaced by a screened interaction, then the self-energy insertion is 
comparable to the G W approximation, which yields generally good quasiparticle band gaps 
(see e.g. [148]), and the exciton line shape is reproduced correctly by the screened electron-hole 
interaction W. (Such a screened kernel will be discussed in section 6.) It has also been shown 
recently [149] that the consistent inclusion of screening beyond EXX in the OEP potential 
reduces the KS eigenvalue gap and brings its value close to the LDA one; the discontinuity 
then yields the correct quasiparticle gap. A proper treatment of correlation (i.e. here: screening) 
is hence crucial for semiconductors, and computational schemes based only on exact exchange 
are not sufficient in this case. In addition, the linearization of the equations as explained in 
section 5.1 is problematic for the gap-opening contribution /(l), because the latter has to 
simulate a discontinuity [32]. This problem can be overcome by using the full, non-linearized 
term as explained in section 6. 



6. Kernels from many-body perturbation theory 

In the previous section the Sham-Schluter equation has been introduced. This equation forms 
a link between many-body perturbation theory (MBPT) and the density-functional framework. 
Since recent results indicate that exploiting this link can be very fruitful, the present section is 
dedicated to the comparison and combination of the two approaches. In fact, for the calculation 
of electronic excitations in solids there is a choice to be made. One can use many-body 
perturbation theory and solve the Bethe-Salpeter equation (BSE), a precise but computationally 
demanding method, or use TDDFT, an in principle computationally more efficient approach, 
but with the limitations of the standard approximations for the exchange-correlation potential 
and kernel (i.e. the TDLDA) outlined above. For an extended review comparing the two 
approaches, see, e.g. [14]. 

The big advantage of TDDFT is that the polarizability is expressed as a variation of the 
density (where the latter is local in space and time); it can be directly obtained from the two- 
point equation (2.35). The BSE uses the more intuitive quasiparticle picture, which makes 
the task to identify efficient approximations easier. However, within MBPT one deals with 
four-point equations. In fact, a key quantity is the four-point reducible polarizability 4 L, which 
can be expressed in terms of the two-particle Green's function G 2 describing the propagation 
of two particles (for absorption, the relevant part describes the propagation of an electron and 
a hole): 



4 L(1, 2, 3, 4) = 4 L (1, 2, 3, 4) - G 2 (l, 2, 3, 4). 



(6.1) 
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As in the previous section, in (6.1) the number '1' represents the set of position, time and 
spin variables (ri, t\,a\). Lo is the disconnected part consisting of two one-particle Green's 
functions G: 

4 L (1, 2, 3, 4) = iG(l, 3)G(4, 2). (6.2) 

The function L and its counterpart L, defined in analogy to x, satisfy a Dyson-like screening 
equation, known as the Bethe-Salpeter equation. When its kernel is approximated to first order 
in the screened Coulomb interaction W it reads 

4 L= 4 L + 4 L ( 4 v- 4 W) 4 L, (63a) 
4 L= 4 L + 4 L ( 4 v- 4 W) 4 L. (63b) 

In (63a) and (63b) we have defined the four-point extension of the Coulomb potential 
Vl, 2, 3, 4) = 5(1, 2)8(3, 4)v(l, 3), whereas 4 W = 8(1, 3)8(2, 4)W(1, 2) is the four-point 
extension of the screened Coulomb potential. Note that the 8 -functions connect different 
indices in the two cases; this is due to the fact that the former stems from variations of 
the Hartree potential, whereas the latter is due to variations of an exchange-like self-energy 
contribution. Because of the way the indices are connected in this second case, the BSE 
cannot be written in a two-point form. Instead, the measurable x is obtained via a two-point 
contraction of 4 L, namely 

P red (l,2) = -L(\, 1,2,2). (6.4) 
From the reduced polarizability P red that is understood to be time (or contour) -ordered, 
the causal response function x can be inferred; the relation between the time-ordered and 
causal response is in fact x(&> > 0) = P Kd (co > 0). Another useful quantity is the 
independent quasiparticle polarizability, Po(l, 2) = — L (l, 1, 2, 2). Note that this is not 
equal to xks = — iGjcs^KS- The expression looks similar, but the Kohn-Sham states and 
eigenvalues in xks are now replaced by their MBPT counterparts, which are, on this level of 
approximation, usually determined in the quasiparticle GW approximation [129] (to be precise, 
most realistic calculations use GW eigenvalues and Kohn-Sham wavefunctions to build G). 

As the two sets of equations (2.35), (3.7), (6.3a) and (63b) have a similar mathematical 
structure, it is natural to try to extract information about the TDDFT exchange-correlation 
kernel through their comparison. Different authors reached very similar expressions for 
the exchange-correlation kernel starting from the BSE but using different approaches [29- 
33,35,146,150]. Tested for real materials, these kernels proved to be successful in reproducing 
the quality of the BSE spectra via a TDDFT formalism. In present implementations their 
computational cost still remains comparable to the cost of solving the BSE but memory 
requirements are significantly reduced, and the expressions can be rewritten leading to 
algorithms with better scaling. Moreover (see section 8), they are a very convenient starting 
point for further approximations. In the following we summarize some of these derivations 
(section 6.2). We will also briefly outline how the combination BSE-TDDFT can be used 
to improve upon current approximations of MBPT, by introducing the TDDFT concept into 
MBPT (section 6.3). Some applications are shown in section 6.4. 



6.1. The exchange-correlation kernel from the Sham-Schliiter equation 

Let us first come back to the non-linearized Sham-Schliiter equation (5.1). As pointed out in 
section 5, it is equivalent to the condition n(r, t) = — iG(r, r' ,t,t + ). In order to obtain an 
exact expression for / xc , Bruneval et al [150] started from this equality and took its functional 
derivative with respect to the density. 8G/8n = — G(8G~ l /8n)G leads to 

, 8G~ l (34) 

iG(13)G(41 + ) — = 8(12) . (6.5) 
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Since the same exact density, and hence the same Hartree potential, should also be obtained 
from the Kohn-Sham potential l>ks = ^ h + v ext + fxc one can write 

G- 1 (12) = G - 1 (12) - 5(12)(uks(1) - v xc (l)) - £(12). (6.6) 

As 8Gq 1 /8n = 0, (6.5) becomes 

P (13)Xks 02) - iG(13)G(41 + )^|^ - P (13)Ac(32) = 5(12). (6.7) 

As in the linearized case (section 5), the exact f xc turns out to consist of two terms 
f xc = f£) + f$ [33, 150]. Those read now as 

fg\\2) = Xks(12) - P _1 (12) (6.8) 



and 



/ x ( c 2) (12) = -iP _1 (l, l')G(l'3)G(41' + )^^. (6.9) 

exactly changes the KS response function into the independent QP one; in particular, 
it therefore solves the band gap problem, f x ^ accounts for the electron-hole interaction. 
Altogether, TDDFT then yields for the irreducible polarizability 

X = Xks + Xks(v + Xks - ^o" 1 + /xc 2) )X = + Po(v + / X ( C 2) )X. (6.10) 

To get an explicit approximation for / xc , one has to choose a starting approximation for the 
self-energy and Green's functions. A simple choice could be to take E, G and Pq as derived 
from a local and adiabatic exchange-correlation potential, e.g. the LDA one. This leads of 
course to the TDLDA and the GWT approach of [151]. A better choice is to start from the 
GW approximation for E, taking W as a screened (e.g. static RPA) Coulomb interaction. For 
the functional derivative, it is then reasonable (i) to neglect the derivative of W as usually 
done in the BSE and (ii) to approximate 8G/8n = -G(8G~ l /8n)G by GP~ l G. The latter 
approximation truncates the chain of derivatives 8H/8H that would appear if one continued 
to calculate all terms of 8G~ l /8n (note that this is equivalent to supposing G be created by a 
local potential). Equation (6.9) then yields 

/4 2 )(34) = j d5d6d7P- 1 (36)G(65)G(5 / 6)W(55 / )G(57)G(75 / )P " 1 (74). (6.11) 

Equation (6.11) is the successful electron-hole exchange-correlation kernel of [29- 
33,35,146,150]. Comparison with (5.4) shows that the two expressions are very similar; 
however, now the formerly bare interaction is screened, which removes the main problems of 
the bare exchange approach. Moreover, QP instead of Kohn-Sham Green's functions appear 
throughout. This reflects the fact that (6.8) has replaced (5.5), which is important since, as 
pointed out above, the linearization of this contribution is in fact problematic. 



6. 2. Comparison of TDDFT and MBPT 

It is instructive to have a look also at the alternative derivations of this kernel. In fact, instead 
of starting from the density it may be more straightforward to start from the observation that 
MBPT and TDDFT should yield the same two-point response function. Since x can in principle 
be obtained via P red from (6.4), one could invert the screening equation and determine f xc from 
f xc = — x~ l — v [152, 153]. The latter relation is of course not of practical interest but 
can give insight on the general features of the kernel, such as its overall frequency behaviour, 
an idea that has been exploited in [154]. 
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Reining et al [29, 30] instead used the equality of the response functions in order to obtain 
an approximate expression for the exchange-correlation kernel of TDDFT by mapping directly 
the matrix elements of the BSE (63a) and (63b) onto the matrix elements of the TDDFT linear- 
response equation (2.35). The derivation starts by rewriting the TDDFT Dyson equation (2.35) 
in a four-point formalism: 

4 X= 4 Xks+ 4 Xks( 4 ^ + Vxc) 4 X, (6.12) 

where the four-point Coulomb interaction 4 u(l,2,3,4) = 8(1,2)8(3, 4)v(l, 3) and the four- 
point kernel 4 f xc (l, 2, 3, 4) = 5(1, 2)8(3, 4)/ xc (l, 3) are defined [30]. 

On the other hand, the BSE has the same structure as (6.12) (cf (63a) and (63b)) but with 
4 Xks replaced by its four-point MBPT counterpart 4 L and the TDDFT kernel replaced by 



4 



K = 8(1, 2)8(3, 4)v(l, 3) - 8(1, 3)8(2, 4)W(1, 2). (6.13) 



In [29, 30], an approximation for both f$ and f$ were derived; since in practice it is 
convenient to use the exact expression (6.8) for the former contribution, only the derivation of 

(2) 

fxc will be discussed in the following. 

We thus replace xks by Po and focus on the second part of the exchange-correlation kernel. 
The first term of the total kernel in both TDDFT and BSE, i.e. the bare Coulomb interaction 

— (2.) 

v, is identical. One might hence try to identify the remaining part of the kernel / xc with the 
screened Coulomb interaction of the BSE kernel (6.13): 

5(1, 2)8(3, 4)f£\l, 3) -8(1, 3)8(2, 4)W(1, 2). (6.14) 

However, as the ^-functions do not contract the same indices expression (6.14) cannot be an 
equality. Different approaches have been proposed to overcome this difficulty and nevertheless 
use the similarity of the equations. 

( i) The assumption that kernels could be similar for a limited number of transitions. 
For this approach, it is useful to write (6.14) on the basis of pairs of Kohn-Sham states 
(transition space) [155]: 

77 TDDFT 77BSE ( r 1 S x 

r {vck)(v'c'k') ^ r {vck){v'c'k'y yv.^J 

with the matrix elements 

F (™c'*') = 2 jd 3 r d 3 r'd>(vkck, r)f£\r, r' , coW(v f k f c'k' , r'), (6.16a) 



F lfkwc>k>) = - f d ' r d 3 r'®(vkv'k\ r)W(r, r' , 0)<t>*(ckc f k f , r'), (6. 



16b) 



where O (ik jk r , r) is the product of a pair of Kohn-Sham wavefunctions; the indices {i, k} stand 
for the band-index and momentum of the Kohn-Sham state. Only the resonant contribution 
is shown (transitions from occupied to empty states) and the momentum transfer is supposed 
to be vanishing since the approach has been derived for optical spectra; both conditions can 
however be easily generalized. As usual in calculations based on the Bethe-Salpeter equation, 
W is moreover supposed to be static. 

Note that (i) in practice, the equality in (6.15) is then imposed for transitions belonging to 
a certain frequency range. The TDDFT spectra will match the BSE spectra only in this region, 
(ii) ^™^J c /£ ) is now a static quantity over that energy range, as a consequence of the static 
approximation for W. (However, since because of (i) a different kernel will be obtained for a 
different group of transitions, there is an overall effective frequency dependence.) 
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Reining et al [29] have formally inverted (6.15). A Fourier transform to reciprocal space 
then leads to 

f£\q -+ 0, G, G') = *Gb*<*; q -+ 0)F^% (v , c , k , ) (<P%}(v'k'c'k'-, q -+ 0). 

VCk v'c'k' 

(6.17) 

The TDDFT exchange-correlation kernel derived in this way has the correct asymptotic 
behaviour (see the discussion in section 8); it stems from the asymptotic behaviour of the 
O, and not from the Coulomb interaction W as one might have suspected. The inversion of 
the matrices O is purely formal; the kernel (6.17) itself was in fact not evaluated in [29] but 
used to derive a TDDFT-like screening equation. A more direct way to arrive at the same 
screening equation along similar lines has been proposed by Sottile et al [30]. One starts by 
writing the TDDFT-like response equation (3.7) in a symmetrized form (the same can be done 
with (2.35)): 

X = Wo - PovPo ~ PoAf^or 1 ^. (6.18) 

It appears that in order to calculate x one on ly needs Pofx?Po and not f$ explicitly. This 
integral contains sums over transitions involving F™^,^ . Enforcing therefore the equality 
in (6.15) yields 

Pof^Po = T 2 , (6.19) 
where T 2 is calculated as 

T 2 (G,G,w) = ^ ^ op _ P QP _,Avckwc*) op _ op _ " ( 6 - 20 > 

k vckv'c'k' £ ck £ vk W 8 c 'k' £ v'k' W 

(here we give again only the expression for the resonant contribution, which is the dominant 
term in absorption spectra). 

This is a result that can be used in practice. The first calculation has been done for bulk 
silicon; the result shown in figure 1 1 obtained by Sottile et al shows almost perfect agreement 
with the result obtained by solving the BSE. 

From a computational point of view, one still has to calculate W; in today's implementation 
moreover the two-particle matrix elements W^xy^') are calculated, which is often the most 
expensive part of a BSE calculation. In principle, the sum over transitions could however 
be performed in a different order and the latter calculation avoided. Moreover by studying 
(6.20), one can derive model kernels, which can drastically decrease the computational cost. 
The performance of the kernel as well as efficient approximations will be illustrated in 
section 8. 

(ii) A perturbative approach leading to the same exchange-correlation kernel was 
proposed by Adragna et al [31 ]. 

Their starting point is again the requirement that the TDDFT and the BSE two-point 
response function should be equal. The polarizabilities are then developed perturbatively 
to the same order. Truncation to first order leads again to (6.19) with (6.20). Marini et al 
then used this kernel for the calculation of optical absorption of insulators, including a bound 
exciton, and loss spectra at non- vanishing momentum transfer [32]. Both cases were in good 
agreement with experiments. It was also shown that the inclusion of the second-order term in 
the perturbative expansion did not significantly affect the calculated spectra. For more details 
see chapters 10 and 20 of [38]. 

(Hi) The perturbative approach was also explored in terms of Feynman diagrams 
[33, 101, 146, 157]. 
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co(eV) 

Figure 11. Imaginary part of the macroscopic dielectric function for Si [30]. Dots: experiment 
[156]; dashed-dotted-dotted line: RPA calculation; dashed-dashed-dotted line: TDLDA 
calculation; dashed-dotted line: GW-RPA calculation; dashed line: BSE calculation; solid line: 
T 2 kernel. 

Stubner, Tokatly and Pankratov [33] derived an integral equation which leads to (6.20) 
when non-locality beyond the first order is neglected. Also the separation of f$ and f$ 
naturally comes about from that derivation. 

Finally we mention the work of von Barth et al [35]: here the kernel is derived from 
a functional approach. Again, this yields in its linearized version as well as the above 

(2) 

approximation for / x y , when GW diagrams are chosen. The advantage of deriving a kernel 
as a double derivative of a functional is that the symmetry of the expression guarantees that 
important conservation laws are fulfilled. This kernel is hence the result of the work of several 
groups, the great majority (although not all) of them being members of a network called 
'nanoquanta' 6 . For convenience, we therefore suggest calling it the 'nanoquanta kernel'. 

6.3. Combining TDDFT and MBPT 

In order to complete the section about links between MBPT and TDDFT, it is useful to recall 
that improved TDDFT response functions and kernels can in turn be inserted into many-body 
calculations in order to go beyond the existing approximations, in particular the RPA-based 
GW approximation. 

The kernel of linear-response TDDFT stems from the variations of the Hartree and 
exchange-correlation potentials with respect to the density. Within MBPT, the variations are 
taken with respect to the one-particle Green's function. However, when one is only interested 
in the self-energy and Green's function (and not in the full four-point response function) in the 
case of MBPT also these variations can be created by a local external potential, for which the 
Runge-Gross theorem holds [16]. As a consequence, in the MBPT scheme it is also possible 
to rely on the fact that the density variations determine the physics of excitations. 

6 http://www.cmt.york.ac.uk/nanoquanta 
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Starting from this observation, Bruneval et al suggested an alternative way to find 
approximations for the self-energy X [150]. This operator, which accounts for all the many- 
body effects on the one-particle Green's function beyond the Hartree term and is a key quantity 
of MBPT, is defined as [129, 130] 

E(l, 2) = i Jd3 d4d5 G(l, 4)r(4, 2; 5)^||i;(3, 1 + ), (6.21) 

where G is the one-particle Green's function, U is a local external potential, V(l) = 
U(l) + l>h(1) is the total classical potential — also including the Hartree contribution t> H , and 
8V(5)/8U(3) = e _1 (5, 3). The irreducible vertex function V is defined as 

8G~ l (l,2) 5£(1,2) 

r(l, 2, 3) = ^— = 5(1, 3)8(2, 3) + — ^— . (6.22) 

8V(3) 8V(3) 

Bruneval et al proposed to replace the chain rule usually employed to obtain 5E/5 V, namely 
(5£/5G)(5G/5V) [130, 158], with the alternative chain rule (5£/5p)(5p/5 V). This step 
is justified by the one-to-one relation between the time-dependent density and the external 
potential or consequently between the density and the classical potential V. Equation (6.22) 
then becomes 

r(l, 2; 3) = 5(1, 3)5(2, 3) + [ d4 ^ (1 ' 2) P(4, 3), (6.23) 

where P = 8p/8V is the irreducible polarizability 

P(l, 2) = -i j d3 d4 G(l, 3)G(4, 1)T(3, 4, 2). (6.24) 

On the other hand, if one multiplies (6.23) with two Green's functions G and integrates, the 
result is 



/ 



P(l, 2) = P (l, 2) + / d3 d4 P (l, 3)/£>(3, 4)P(4, 2) (6.25) 
with (6.9). 

This constitutes therefore an alternative derivation of the exact kernel. Moreover, once 
P is known (6.23) allows one to calculate the three-point vertex, and hence an improved self- 
energy, without solving a four-point integral equation. It has been shown [150] that, besides 
the better description of the test-charge-test-charge screened W, a major contribution to the 
correction beyond GW stems from the induced exchange-correlation potential that acts on an 
additional particle or hole on top of the induced Hartree potential; only the latter is contained 
in GW (moreover, it is approximated in RPA). 



6.4. Applications 

When considering practical calculations it should be noted that in general (i) as mentioned 

(2) 

above, only / x v c is explicitly used, whereas the GW correction is included by using Pq, (ii) most 
calculations are restricted to the resonant contribution and (iii) the diagonal contribution 
F(ldc)(vck) * s extrac t e d from /x? and included in Po. Strictly speaking, this contribution should 
be zero in a solid, but for a finite &-point sampling it is non- vanishing. These three points are 
important for the quality of the final results. 

An example of the application of the BSE-derived kernel (6.20) to real systems can be 
found in figure 12 [32] , where the optical absorption spectrum of LiF is shown. One can observe 
that, as anticipated, the agreement between the BSE calculation and the TDDFT calculation 
using the BSE-derived kernel is excellent. The importance of the matrix character of f xc 
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Figure 12. Optical absorption spectra calculated within the BSE (dots), and TDDFT (solid line) [32] 
compared with experiments (circles) [88]. The results obtained using a scalar f xc (dashed line) are 
shown to stress the importance of the G, G' / terms in / xc . 
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Figure 13. Optical absorption spectra of polyacethylene compared with experiment (vertical dashed 
line) [161, 162]. The BSE calculation (dots) and the TDDFT calculation using the many-body f xc 
(solid line) are in excellent agreement with the experiment. GW-RPA (dashed line) and TDLDA 
calculations (dotted-dashed line) are also shown. 



(i.e. the importance of G, G' ^ terms) in wide-gap insulators is demonstrated by looking at 
the curve obtained with a scalar (G = G' = 0) f xc . Although the kernel is strongly frequency 
dependent, only the main peak is reproduced correctly, while some unphysical regions of 
negative absorption appear at higher energies. In contrast, when dealing with systems with 
continuum exciton effects (Si, diamond, etc) the use of the head of f xc is often sufficient to 
recover the BSE spectrum. This is related to the different degree of inhomogeneity of the 
induced density in simple semiconductors and wide-gap insulators. 

The BSE-derived kernel (6.20) has also been successfully applied to the study of low- 
dimensional systems. An example is conjugated polymers. We have already pointed out 
in section 3.2 that simple adiabatic local-density and gradient-corrected functionals fail in 
describing the dielectric response of long molecular chains [24]. Current density-functional 
approaches can partially solve this problem [159, 160]. The non-local frequency-dependent 
BSE-derived kernel restores the good agreement between the calculated and the measured 
polarizabilities. This can be observed in figure 13, where the calculated optical absorption 
spectrum of polyacethylene is compared with experimental data [161, 162]. 
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In brief, the BSE-derived f xc kernel is able to reproduce the optical and energy-loss spectra 
of a large class of materials including semiconductors and large band gap insulators [30, 32] 
as well as systems of low dimensionality, including surfaces [161]. Although calculations 
today are still relatively cumbersome, improved algorithms and/or efficient approximations 
are expected to make this approach clearly competitive with respect to a BSE calculation, 
when one is interested in an efficient determination of spectra. 

7. Current-density functionals 

Finally, we would like to briefly outline an alternative way to approach the question of 
going beyond the ALDA. In fact, many problems of the latter, including its dramatic failure 
for optical absorption spectra of insulators, are related to the fact that it is derived only 
from the local density. As a consequence, the effective potential inside a bulk solid is 
insensitive to the accumulation of surface charges, which may have a large influence on the 
measured material properties in actual experiments. In particular, they give rise to a long- 
range exchange-correlation potential that counteracts an applied external field, thus reducing 
the observable polarizability. As the density inside the bulk does not change in this case, 
such a counteracting term cannot be reproduced by a purely local approximation. Instead, 
the exchange-correlation potential must depend in a non-local manner on the global density 
distribution. This is difficult to achieve with explicit density functionals, which are, like the 
ALDA, typically derived from the homogeneous electron gas. Implicit density functionals 
constructed from the Kohn-Sham wavefunctions do not suffer from the same limitation, but 
the large computational cost makes it difficult to incorporate high-order correlation terms in 
practice. An alternative approach to include non-locality is given by time-dependent current- 
density-functional theory, because, in contrast to the density n(r, t), the current density j(r, t) 
can be used as a local indicator of global changes in the system. For example, its behaviour 
correctly reflects the charge transport through the bulk between two solid surfaces upon 
polarization [163]. 

Current-density-functional theory was originally developed as a generalization of TDDFT 
for time-dependent magnetic fields [56, 164]. In this case, in addition to the scalar potential 
fext(/\ t), the Hamiltonian 



contains an external vector potential A ext (r, t). The conjugate variable that couples to 
the latter is the current density, defined in an analogy to (2.5) but with — iV + A ext (r, t) 
instead of the canonical momentum operator — iV. Within this framework the auxiliary 
Kohn-Sham system is chosen in such a way that it yields not only the same density but 
also the same current density as the original interacting electron system. As a necessary 
prerequisite for achieving non-interacting f-representability of both quantities, the effective 
vector potential A e ff (r, t) = A ext (r, t) + A xc (r, t) in the Kohn-Sham equations must then 
feature a non-vanishing exchange-correlation contribution, because changes in the scalar 
potential alone have no effect on the transverse component of the current [165]. The 
scalar potential can actually be eliminated at all times through the gauge transformation 
A ext (>, t) -> A ext (r, t) + VA(r, t) and v QXt (r, t) -> v QXt (r, t) - 3A(r, t)/dt, if the gauge 
function satisfies 3A(r, t)/dt = v QXt (r, t) with the initial condition A(r, to) = 0. Besides, 





(7.1) 
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the evolution of the density is completely determined through the static equilibrium value 
n(r,to) = riQsir) and the continuity equation (2.8). Therefore, current-density-functional 
theory can be formulated solely in terms of vector potentials and with the current density as 
the only independent variable. An elegant proof of their one-to-one correspondence, which 
includes the Runge-Gross theorem as a special case, was given by Vignale [166]; the vector 
potential is determined uniquely up to a gauge transformation that does not alter the initial 
state. 

For practical purposes the exchange-correlation vector potential must be specified as 
a functional of the current density. Fortunately, it turned out that a consistent gradient 
expansion around the limit of the homogeneous electron gas is possible in this approach. 
By investigating weakly inhomogeneous systems, Vignale and Kohn [167] thus derived an 
approximation for the induced field that includes dynamic exchange and correlation effects 
and depends only on the local current density. The resulting expression can be interpreted 
in terms of viscoelastic stresses in the electron fluid [168, 169]. So far the approximation of 
Vignale and Kohn has mainly been applied to plasmon line widths in semiconductor quantum 
wells within the one-band effective-mass model [170-172] as well as to one-dimensional 
polymers and other extended macromolecules, where it leads to significant improvements in 
the polarizabilities over a treatment within the ALDA [160, 173]. A simplified version has also 
been used in the context of solid-state physics to calculate optical absorption spectra [159]. 
The results reproduce the essential spectral features of several common semiconductors and 
insulators, including the excitonic E\ resonance in silicon, and are in much better agreement 
with experimental data than the ALDA. However, the absorption threshold is too low, and it 
remains unclear whether the residual discrepancies are due to the Vignale-Kohn approximation 
itself or the additional simplifications of the implementation. It should be noted that the 
simplified formula applied in [159] is very similar to the long-range approximation [29, 34] 
that will be discussed in the next section. This explains both its success in describing 
the line shape and its failure to open the band gap without the help of an explicit energy 
shift (a task that is carried out by the GW corrections in the case of the long-range kernel 
of section 8). 

8. Simple models 

In section 6 we introduced a class of parameter-free kernels that are particularly successful 
for the description of optical absorption in solids — the kernels derived from the BSE. 
Even though they have potentially reduced computational effort with respect to the BSE, 
calculations using these kernels are however significantly more involved than those within 
the RPA or the TDLDA. Therefore, the question of finding simple and efficient, but also 
reliable, models for f xc is still open. In this quest, there are several important lessons to 
be learnt from the BSE-derived kernels. In the following we discuss some model kernels, 
inspired by the BSE, that combine simplicity with a reasonably good description of response 
properties. 

8.1. Long-range exchange -correlation kernels 

One of the most striking characteristics of the TDLDA kernel is that it is static and local in 
space. Therefore, one can expect that the inclusion of either dynamical (memory) effects or 
long-range non-local terms [21, 146] or both improves, in principle, the results yielded by the 
simple TDLDA. In this section, we introduce some model kernels, obtained by approximating 
the BSE-derived kernel of (6.9) and (6.19), which accounts for such further terms. 
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Figure 14. Material dependence of a statlc with respect to the inverse of the dielectric constant [34]. 



The exchange-correlation kernel (6. 9), (6. 19) contains a long-range contribution (LRC) 
of the form l/q 2 . This LRC is instead completely absent within the TDLDA, as the TDLDA 
kernel tends to a constant in the limit q —> 0. The simplest model that exhibits the correct 
LRC has the form 

a static 

/xc atiC (?) = 2-' ^ 

Q 

where a statlc is a material-dependent parameter. The use of this particular form can be motivated 
starting from the BSE, following the lines of section 6 [29]. In fact, the function <&c=o(vkck; q) 
in (6.17) tends to zero as O ~ q for q — >► 0. Since in this limit F* s f., . , 1f . behaves like a 

v 7 11 (vck)(v'c'k) 

constant, this implies that f xc (q,G = 0, G' = 0) tends to zero as l/q 2 in the optical limit. 
Moreover, as can be seen in figure 14 the LRC to the exchange-correlation kernel is inversely 
proportional to the macroscopic dielectric constant €oq [34] . 

Note that there is also a positive long-range contribution to the exchange-correlation 
kernel stemming from the QP shift of eigenvalues (as predicted in [174]), which competes 
with the negative one resulting from the BSE, i.e. the electron-hole interaction. This 
contribution is contained in (6.8). In [29, 34] it was shown that the LRC alone is sufficient to 
reproduce the contribution of the electron-hole interaction / x v c to the optical spectra of simple 
semiconductors. The same is not true concerning the self-energy contribution where a 
LRC approximation does not work. For this reason in the following discussion we will focus 

(2) 

on models for the term /x C . 

Figure 15 shows the optical absorption spectra of GaAs [34]. As previously discussed 
in section 3.3, the TDLDA result is close to the RPA curve, and both show well-known 
discrepancies with experiment: the peak positions are redshifted and the intensity of the 
first main structure is strongly underestimated. The dashed curve (GW-RPA) is obtained by 
replacing the KS eigenvalues with GW quasiparticle energies in the RPA calculation of €. This 
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Figure 15. Imaginary part of the macroscopic dielectric function of GaAs [34]. Dots: experiment 
[93]; dotted-dashed curve: TDLDA; dashed curve: GW-RPA; solid line: TDDFT-LRC. 




24 



Energy (eV) 

Figure 16. Imaginary part of the macroscopic dielectric function for LiF [87]. Dots: experiment 
[88] ; dotted curve: BSE; dotted-dashed curve: TDLDA; dashed curve: TDDFT using the static LRC 
model kernel; solid line: TDDFT using the dynamical LRC model kernel. In the inset: GW-RPA. 

corresponds to applying the Dyson equation (2.35) for the first part of the exchange-correlation 
kernel f£\ The resulting spectrum is blueshifted and, moreover, the lineshape has not been 
corrected. Finally, the curve representing a TDDFT calculation starting from GW quasiparticle 
energies to calculate the independent-particle response and using the LRC kernel (8.1) gives 
an excellent fit to the experiment. This curve was obtained using a statlc = 0.2. 

The simple static LRC model, together with the linear dependence of a statlc on l/e^, 
allows one to predict the absorption spectrum from the knowledge of the experimental dielectric 
constant of the material in question. However, the model / x s <! atlc of (8.1) has clear limits that 
become more and more evident as the band gap increases. 

For example, in diamond [34, 154] it is not possible to describe with good precision both 
the first shoulder and the main peak in the spectrum using a single parameter a statlc . The 
problem gets even more serious when bound exciton peaks appear in the spectrum as, e.g. 
for LiF [87] (see figure 16) or solid argon [175]. In fact, the only possible action of the 
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LRC exchange-correlation kernel is to redistribute oscillator strength. In contrast to the full 
BSE kernel, this remains also true in the cases where poles should appear in the bandgap, 
such as in LiF where bound excitons occur in the experimental spectrum (see figure 16). 
Furthermore, approximations as outlined above are based on a limited range of transitions. 
When the spectral range gets larger, e.g. when the response beyond the absorption region is 
considered, as a consequence of the change in the limited range of transitions f xc has to change 
also. In the case of the model kernel (8.1) this introduces an effective ^-dependence of the 
parameter a statlc . In order to describe the plasmon of silicon one could for example still use 
(8.1), but with an a statlc that is an order of magnitude larger than the one that yields a good 
optical spectrum [34]. This will be further discussed below, after a brief summary of other 
approaches leading to a LRC. 

The simple static LRC model (8.1) does not in fact represent the first attempt to account 
for long-range effects in semiconductors. In 1994 Godby and Sham [176] pointed out that 
long-range density variations give rise to an effective exchange-correlation field. It was 
then proved [23] that the origin of this exchange-correlation field lies in the macroscopic 
polarization. Alternatively, one could avoid the use of macroscopic polarization by introducing 
a scissor-operator quasiparticle correction to the Kohn-Sham gap. The LRC stemming from 
this discussion has hence the task of increasing the gap (like f£?) and has, consequently, a 
positive sign. 

Aulbur, Jonsson and Wilkins [177] studied the problem of determining the static dielectric 
constants. They related an effective exchange-correlation field to the difference between 
the true and the Kohn-Sham static susceptibilities and, by using calculated Kohn-Sham and 
measured (i.e. 'true') values, fixed the prefactor of a LRC to the kernel for a series of materials. 
It resulted in a contribution Af xc = y/q 2 , where y is of the order of 0.25 for several small and 
medium-gap semiconductors. According to this work, the LRC correction should account for 
both /x<P and which is the reason why y turns out to be small and positive. However, 
this model cannot be extended to the description of spectra, as the quasiparticle corrections 
simulated by f$ are too complex to be written as a simple LRC model. 

Another approach to the description of the dynamical susceptibility was proposed by 
Boeij et al [159] who obtained a polarization-dependent functional derived from current 
density-functional theory [178] (see section 7). This functional involves two parameters: one 
(material-dependent) accounting for a positive shift of transition energies (in fact, a scissor 
operator) and a second (constant) one, chosen to be 0.4, that multiplies a tensor Y containing 
the polarization effects. The tensor Fis in principle frequency-dependent, but its static value 
is used. This approach yields identical results to the static LRC model described here, if one 
identifies 0.4F = — a statlc and if the scissors are replaced by the quasiparticle corrections. 

In consideration of the limits of the static LRC model, it is natural trying to guess how 
to further improve the kernel. One can either work along the lines of [29] and add a more 
complicated spatial behaviour or one can keep the simple l/q 2 form of the LRC, but introduce 
a frequency-dependent a [154]. Clearly, this latter choice allows for each structure of the 
spectrum to have its own effective correction. Botti et al [87] proposed a frequency-dependent 
LRC kernel of the form 

f^(g,a>) = ~(a + Pa?). (8.2) 

q 

This choice was guided by recent calculations [154] for bulk silicon and diamond, which 
yielded the frequency dependence of the LRC term of the exchange-correlation kernel from 
the inversion of the BSE. Figure 17 shows, for f xc of silicon, that the dynamical model / x d c yn is 
indeed a good approximation in a large energy range, even including the plasmon region (see 
inset of figure 17). 
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Figure 17. The real part of the LRC component of the exchange-correlation kernel needed to 
reproduce the BSE optical spectrum of silicon. The inset shows the same function in an energy 
region close to the plasma frequency. Dashed lines: calculation from [154]; solid line: dynamical 
LRC model. 



Starting from (6.20) it is possible to prove that the two parameters a and of (8.2) can 
be related to physical quantities, namely the dielectric constant and the plasma frequency. For 
simple semiconductors, this model yields the same result as the static LRC model [87]. In 
the case of wide-gap insulators, such as diamond or cadmium selenide, or for the EELS of 
silicon, the dynamical model rectifies the shortcomings of the static LRC kernel [87]. The 
improvement is significant even in the presence of bound excitons. We consider again as an 
example LiF (see figure 16). Setting a statlc to the value 2.0 gives a reasonable compromise 
(dashed line), enhancing slightly the low energy structures without provoking the collapse of 
the spectrum. The worst disagreement concerns the absence of the large excitonic peak at about 
12.5 eV. On the other hand, the dynamical model (solid line) is able to describe the strong bound 
exciton peak. A better agreement can only be found using the BSE approach [89, 91, 92, 179] 
(dotted line) or the full kernel derived by the BSE (see figure 12), both of which involve much 
larger computational effort than the application of the simple dynamical LRC model. Note that 
the objective of these model kernels is not to provide an ab initio approximation to the exact 
TDDFT kernel, but to offer a numerically efficient framework for the calculation of response 
properties of complex solids and nanostructures. A key feature is that the model kernels depend 
on a reduced number of parameters and that these can be related in a straightforward manner 
to known physical quantities. The drawback of these models is, of course, the reduced range 
of validity of the underlying approximation, which implies that one should carefully check the 
applicability of a model to a specific system. 



8.2. Contact exciton 

The contact exciton model was successfully used in the 1970s for the description of continuum 
exciton effects in a wide range of systems [180-183]. Model calculations could also show that 
the contact exciton is able to produce one bound state [183, 184]. This simple model gives rise 
to the equation 

€ M (co) - 1 = RPA -. (8.3) 
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Figure 18. Imaginary part of the macroscopic dielectric function for bulk silicon. Dotted 
curve: BSE; Dotted-dashed curve: TDDFT with static LRC (a statlc = 0.2); continuous curve: 
TDDFT contact kernel (A/2 = 15); dots: experiment [156]. 

Inspection reveals that this is the LRC kernel approximation derived in the previous section 
[154, 175], provided that f xc = —gAn/q 1 . Moreover, in [180] the parameter g is dependent 
on co 2 , similarly to the parameter /3 of the dynamical LRC model in section 8.1. 

A different exchange-correlation kernel can be derived in the same spirit by replacing the 
screened electron-hole interaction in the BSE by a local contact potential [175]. When W in 
(6.14) is approximated by 8 (12) A, the equality is possible in (6.15) and one obtains 

/ x c c ontaCt (n,r 3 ) = - i«(n -r 3 )A. (8.4) 

This is obviously an ultra- short-range kernel. 

In figure 18 we compare optical absorption calculations performed within the different 
approximations with the experimental optical spectrum of silicon. The curve obtained for 
the ultra-short-ranged / x c c ontact using A/2 = 15 (continuous curve) turns out to be in very 
good agreement with the result of a full BSE calculation (dotted curve) and the long-range 
kernel f^ atlc = —a/q 2 with a = 0.2 (dotted-dashed curve). All calculations are also in good 
agreement with experiment (circles). 

In order to understand how a short-range and a long-range f xc can yield similar spectra 
for simple semiconductors, one has to realize that we are talking about effective exchange- 
correlation kernels and not the real one. These kernels are chosen to reproduce a certain 
number of properties. It is clear that more than one kernel (even with very different analytical 
forms and physical interpretation) can lead to the same optical spectrum. The contact exciton 
is a very good illustration of this. 

To a certain extent, a short-ranged but strong electron-hole Coulomb interaction can 
simulate the effect of the true, screened long-ranged one: continuum exciton effects can be 
well reproduced, and even one bound exciton can be created. However, in the latter case the 
continuum is not well described, and no higher-order peaks of the series are obtained. 

The choice of which model to use should finally depend both on computational 
convenience and on the possibility of determining the parameters of the model without fitting 
to experiment. The stronger the requirements on the precision of the results, the closer the 
chosen model should resemble the exact f xc and, for example, have a long-range component 
when the system is a semiconductor or insulator. We remark that all these methods discussed 
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above are the result of quite recent investigations, so it is very reasonable to expect further 
developments in the near future. 

9. Conclusions 

After the pioneering work of Zangwill and Soven [15], time-dependent density-functional 
theory (TDDFT) found a rigorous foundation in the work by Runge and Gross [16]. The 
time-dependent Kohn-Sham equations are obtained as a generalization of the static case and, 
from these, the response functions describing the neutral excitations of the system. The main 
ingredient is the time-dependent exchange-correlation potential v xc [n](r, t), which depends 
on the density at all points in space and at all past times. 

Staying within linear-response theory, one needs to also know, in addition to the static 
exchange-correlation (xc) potential f X cl/*GS the so-called xc kernel f xc [n](r,r f , t,t f ) = 
8v xc [n](r, t)/8n(r f , t r ). If these two quantities are known, TDDFT is an exact theory, yielding 
the exact linear response. The problem is, therefore, how to generate suitable approximations 
for the potential and the kernel. For finite systems, it is sometimes crucial to have a good v xc , 
whereas in extended systems in general the main problem is to find a good approximation for 
the kernel f xc . 

Neutral excitations can also be calculated in the framework of many-body perturbation 
theory, via the solution of the Bethe-Salpeter equation (BSE). However, in that case one 
has to deal with two-particle Green's functions, which are four-point quantities, whereas 
TDDFT is based on two-point linear response functions. Therefore, TDDFT promises to 
be computationally much more efficient, which motivates the search for good approximations. 

The simplest approximation is the adiabatic LDA (TDLDA), which often yields good 
results for finite systems. Its main shortcoming is to miss the long-range part, proportional 
to l/\r — r'|, which may be important in extended systems. Nevertheless, the TDLDA also 
describes well some properties of extended systems, in particular, the plasmon structures in the 
loss function or the dynamical structure factor. In these cases, indeed, the long-range part of the 
kernel is not important, either because of cancellations or because of non- vanishing momentum 
transfer. However, TDLDA does not describe well the optical properties of extended systems, 
where the long-range part is essential. 

Therefore, one of the main challenges of recent years has been to find suitable 
approximations to the long-range part of f xc . Promising solutions came from a comparison of 
TDDFT and the many-body perturbation theory (MBPT) equations for polarization. Different 
groups have suggested different approaches; all approaches have led to one and the same 
kernel, the 'nanoquanta' kernel, that is linear in the screened Coulomb interaction and based 
on quasiparticle ingredients. This kernel allows one to obtain a good description of the optical 
properties of extended systems, including bound excitons. It also describes well the spectra 
for non-vanishing momentum transfer, and it has been demonstrated that it yields very good 
results for low-dimensional systems, including clusters, wires and surfaces. 

At present, although the TDDFT linear-response equation is relatively quick to solve 
because of its two-point nature, calculations involving the 'nanoquanta' kernel are still 
computer-time intensive. In particular, for the calculation of the kernel itself current implemen- 
tations use ingredients of the BSE approach, so that the performance of the two methods is com- 
parable. However, the current optimization of algorithms and computer codes should lead to 
considerable speedup of the TTDFT calculations. Moreover it has been shown that the new ker- 
nel is a convenient starting point for additional approximations, such as the LRC a/q 2 method 
described in section 8.1. Hence, we believe that TDDFT using the many-body kernel will be 
one of the main methods in the future for calculating the optical properties of extended systems. 



404 



S Botti et al 



Acknowledgments 

This work was funded in part by the EU's 6th Framework Programme through the Nanoquanta 
Network of Excellence (NMP-4-CT-2004-500198). The authors acknowledge discussions and 
collaborations with M Gatti, S Kurth, A Marini, V Olevano, O Pulci, A Rubio, M Scheffler, 
F Sottile, K Tatarczyk, I Tokatly and U von Barth. They thank I G Gurtubay, A Marini, 
G Mulas and D Varsano for providing their figures included in this paper. They also thank 
M Marques for the critical reading of the manuscript. RDS acknowledges support from PRIN 
2005 ofMIUR. 



References 

[1] Kohn W 1999 Rev. Mod. Phys. 71 1253-66 

[2] Hohenberg P and Kohn W 1964 Phys. Rev. 136 B 864-71 

[3] Kohn W and Sham L J 1965 Phys. Rev. 140 Al 133-8 

[4] Perdew J P, Burke K and Ernzerhof M 1996 Phys. Rev. Lett. 11 3865-8 

[5] Perdew J P, Burke K and Ernzerhof M 1997 Phys. Rev. Lett. 78 1396(E) 

[6] Jones R O and Gunnarsson O 1989 Rev. Mod. Phys. 61 689-746 

[7] Baroni S and Resta R 1986 Phys. Rev. B 33 7017-21 

[8] Kotliar G, Savrasov S Y, Haule K, Oudovenko V S, Parcollet O and Marianetti C A 2006 Rev. Mod. Phys. 
78 865-951 

[9] Sim F, St-Amand A, Papai I and Salahub D R 1992 J. Am. Chem. Soc. 114 4391-400 
[10] Godby R W, Schluter M and Sham L J 1987 Phys. Rev. B 36 6497-500 

[11] Fetter A L and Walecka J D 1971 Quantum Theory of Many-Particle Systems (New York: McGraw-Hill) 
[12] Sham L J and Rice T M 1966 Phys. Rev. 144 708-14 
[13] Hanke W 1978 Adv. Phys. 21 287-341 

[14] Onida G, Reining L and Rubio A 2002 Rev. Mod. Phys. 74 601-59 and references therein 

[15] Zangwill A and Soven P 1980 Phys. Rev. Lett. 45 204-7 

[16] Runge E and Gross E K U 1984 Phys. Rev. Lett. 52 997-1000 

[17] Gross E K U and Kohn W 1985 Phys. Rev. Lett. 55 2850-2 

[18] Gross E K U and Kohn W 1986 Phys. Rev. Lett. 57 923(E) 

[19] van Leeuwen R 2001 Int. J. Mod. Phys. B 15 1969-2023 

[20] Petersilka M, Gossmann U J and Gross E K U 1996 Phys. Rev. Lett. 76 1212-5 

[21] Marques M A L and Gross E K U 2004 Annu. Rev. Phys. Chem. 55 427-55 

[22] Petersilka M, Gross E K U and Burke K 2000 Int. J. Quantum Chem. 80 534-54 

[23] Gonze X, Ghosez P, and Godby R W 1995 Phys. Rev. Lett. 74 4035-8 

[24] van Gisbergen S J A, Kootstra F, Schipper P R T, Gritsenko O V, Snijders J G and Baerends E J 1999 

Phys. Rev. Lett. 83 694-7 
[25] Maddocks N E, Godby R W and Needs R J 1994 Europhys. Lett. 21 681-6 

[26] Gurtubay I G, Ku W, Pitarke J M, Eguiluz A G, Larson B C, Tischler J and Zschack P 2005 Phys. Rev. B 
72 125117 

[27] Waidmann S, Knupfer M, Arnold B, Fink J, Fleszar A and Hanke W 2000 Phys. Rev. B 61 10149-53 

[28] Weissker H C, Serrano J, Huotari S, Bruneval F, Sottile F, Monaco G, Krisch M, Olevano V and Reining L 

2006 Phys. Rev. Lett. 91 237602 
[29] Reining L, Olevano V, Rubio A and Onida G 2002 Phys. Rev. Lett. 88 066404 
[30] Sottile F, Olevano V and Reining L 2003 Phys. Rev. Lett. 91 056402 
[31] Adragna G, Del Sole R and Marini A 2003 Phys. Rev. B 68 165108 

Del Sole R, Pulci O, Olevano V and Marini A 2005 Phys. Status Solidi b 242 2729-36 
[32] Marini A, Del Sole R and Rubio A 2003 Phys. Rev. Lett. 91 256402 
[33] Stubner R, Tokatly I V and Pankratov O 2004 Phys. Rev. B 70 2451 19 

[34] Botti S, Sottile F, Vast N, Olevano V, Weissker H C, Reining L, Onida G, Rubio A, Del Sole R and Godby R W 

2004 Phys. Rev. B 69 155112 
[35] von Barth U, Dahlen N E, van Leeuwen R and Stefanucci G 2005 Phys. Rev. B 72 235109 
[36] Ankudinov A L, Nesvizhskii A I and Rehr J J 2003 Phys. Rev. B 67 1 15120 

[37] Kootstra F, de Boeij P L, van Leeuwen R and Snijders J G 2002 Reviews in Modern Quantum Chemistry: A 
Celebration of the Contributions of Robert Parr ed K D Sen (Singapore: World Scientific) 1155-85 



Time-dependent density-functional theory for extended systems 



405 



[38] Marques M A L, Ullrich C A, Nogueira F, Rubio A, Burke K and Gross E K U (ed) 2006 Time-Dependent 

Density Functional Theory (Lecture Notes in Physics) (Heidelberg: Springer) 
[39] Nunes R W and Vanderbilt D 1994 Phys. Rev. Lett. 73 712-5 
[40] Maitra N T and Burke K 2001 Phys. Rev. A 63 042501 
[41] Maitra N T and Burke K 2001 Phys. Rev. A 64 039901(E) 
[42] van Leeuwen R 1999 Phys. Rev. Lett. 82 3863-6 
[43] Levy M 1982 Phys. Rev. A 26 1200-8 
[44] Lieb E H 1983 Int. J. Quantum Chem. 24 243-77 
[45] van Leeuwen R 1998 Phys. Rev. Lett. 80 1280-3 

[46] Gross E K U, Ullrich C A and Grossman U J 1995 Density Functional Theory (NATO ASI Series B, vol 337) 

ed E K U Gross and R M Dreizler (New York: Plenum) pp 149-71 
[47] Gross E K U, Dobson J F and Petersilka M 1996 Density Functional Theory (Topics in Current Chemistry vol 

181), ed R F Nalewajski (Heidelberg: Springer) p 81 
[48] Keldysh L V 1965 Sov. Phys.—JETP 20 1018-26 

[49] Ullrich C A, Gossmann U and Gross E K U 1995 Phys. Rev. Lett. 74 872-5 
[50] Gorling A 1997 Phys. Rev. A 55 2630-9 
[51] http://theory.polytechnique.fr/codes/dp 
[52] http://www.fisica.uniroma2.it/ self 
[53] http : // exciting . sourcef orge .net 

[54] Bertsch G F, Iwata J I, Rubio A and Yabana K 2000 Phys. Rev. B 62 7998-8002 

[55] Castro A, Marques M A L and Rubio A 2004 /. Chem. Phys. 121 3425-33 

[56] Ghosh S K and Dhara A K 1988 Phys. Rev. A 38 1 149-58 

[57] Kootstra F, de Boeji P L and Snijders J G 2000 J. Chem. Phys. 112 6517-31 

[58] Walker B, Saitta A M, Gebauer R and Baroni S 2006 Phys. Rev. Lett. 96 1 13001 

[59] Gunnarsson O and Lundqvist B I 1976 Phys. Rev. B 13 4274-98 

[60] Gunnarsson O and Lundqvist B I 1977 Phys. Rev. B 15 6006(E) 

[61] Perdew J P and Levy M 1985 Phys. Rev. B 31 6264-72 

[62] Ehrenreich H and Cohen M 1959 Phys. Rev. 115 786-90 

[63] Adler S L 1962 Phys. Rev. 126 413-20 

[64] Wiser N 1963 Phys. Rev. 129 62-9 

[65] Del Sole R and Fiorino E 1984 Phys. Rev. B 29 4631-45 

[66] Sottile F, Bruneval F, Marinopoulos A G, Dash L K, Botti S, Olevano V, Vast N, Rubio A and Reining L 2005 

Int. J. Quantum Chem. 102 684-701 
[67] Marques M A L, Castro A and Rubio A 2001 /. Chem. Phys. 115 3006-14 

[68] Perdew J P and Kurth S 2003 A Primer in Density Functional Theory (Lecture Notes in Physics vol 706) 

ed C Fiolhais, F Nogueira and M Marques (Berlin: Springer) pp 1-55 
[69] Hybertsen M S and Louie S G 1986 Phys. Rev. B 34 5390-413 
[70] Lindhard J 1954 K. Dan. Vidensk. Selsk. Mat.-Fys. Medd. 28 1-57 
[71] Perdew J P and Burke K 1996 Int. J. Quantum Chem. 57 309-19 
[72] Perdew J P, Kurth S, Zupan A and Blaha P 1999 Phys. Rev. Lett. 82 2544-7 
[73] Perdew J P, Kurth S, Zupan A and Blaha P 1999 Phys. Rev. Lett. 82 5179(E) 
[74] Castro A, Marques M A L, Alonso J A and Rubio A 2004 /. Comput. Theor. Nanosci. 1 231-55 
[75] Gurtubay I G, Pitarke J M, Ku W, Eguiluz A G, Larson B C, Tischler J, Zschack P and Finkelstein K D 2004 

Phys. Revolts 201201(R) 
[76] Ando T 1977 Z. Phys. B 26 263-334 

[77] Yabana K and Bertsch G F 1999 Int. J. Quantum Chem. 75 55-66 

[78] Onida G, Schmidt W G, Pulci O, Palummo M, Marini A, Hogan C and Del Sole R 2001 Phys. Status Solidi a 
188-1242 1233-42 

[79] Marques M A L, Lopez X, Varsano D, Castro A and Rubio A 2003 Phys. Rev. Lett. 90 258101 
[80] Marques M A L and Botti S 2005 /. Chem. Phys. 123 014310 

[81] Rappoport D and Furchte F 2006 Time-Dependent Density Functional Theory (Lecture Notes in Physics vol 

706) ed M A L Marques et al (Berlin: Springer) pp 337-54 (chapter 23) 
[82] Wang C R C, Pollack S, Cameron D and Kappes M M 1990 J. Chem. Phys. 93 3787-801 
[83] Malloci G, Mulas G and Joblin C 2004 Astron. Astrophys. 426 105-1 17 
[84] Joblin C, Leger A and Martin P 1992 Astrophys. J. Lett. 393 L79-81 

[85] Champagne B, Perpete E A, van Gisbergen S J A, Baerends E J, Snijders J G, Soubra-Ghaoui C, Robins K A 

and Kirtman B 1998 /. Chem. Phys. 109 10489-98 
[86] Gritsenko O V and Baerends E J 2001 Phys. Rev. A 64 042506 



406 



S Botti et al 



[87] Botti S, Fourreau A, Nguyen F, Renault Y O, Sottile F and Reining L 2005 Phys. Rev. B 72 125203 

[88] Roessler D M and Walker W C 1967 /. Opt. Soc. Am. 57 835-6 

[89] Albrecht S, Reining L, Del Sole R and Onida G 1998 Phys. Rev. Lett. 80 4510-3 

[90] Benedict L X, Shirley E L and Bohn R B 1998 Phys. Rev. Lett. 80 4514-7 

[91] Rohlfing M and Louie S G 1998 Phys. Rev. Lett. 81 2312-15 

[92] Rohlfing M and Louie S G 2000 Phys. Rev. B 62 4927-44 

[93] Lautenschlager P, Garriga M, Vina L and Cardona M 1987 Phys. Rev. B 36 4821-30 

[94] Stiebling J 1978 Z. Phys. B 31 355-7 

[95] Olevano V and Reining L 2001 Phys. Rev. Lett. 86 5962-5 

[96] Marinopoulos A, Reining L, Olevano V, Rubio A, Pichler T, Liu X, Knupfer M and Fink J 2002 Phys. Rev. Lett. 
89 076402 

[97] Vast N, Reining L, Olevano V, Schattschneider P and Jouffrey B 2002 Phys. Rev. Lett. 88 037601 

[98] Dash L K, Vast N, Baranek P, Cheynet M C and Reining L 2004 Phys. Rev. B 70 245116 

[99] van Leeuwen R and Baerends E J 1994 Phys. Rev. A 49 2421-31 

[100] Dobson J F, Biinner M J and Gross E K U 1997 Phys. Rev. Lett. 79 1905-8 

[101] Tokatly I V 2006 Time-Dependent Density Functional Theory (Lecures Notes in Physics vol 706) 

ed M A L Marques et al (Berlin: Springer) pp 123-36 (chapter 8) 

[102] Maitra N T 2005 /. Chem. Phys. 122 234104-9 

[103] Maitra N T, Zhang F, Cave R J and Burke K 2004 /. Chem. Phys. 120 5932-7 

[104] Cave R J, Zhang F, Maitra N T and Burke K 2004 Chem. Phys. Lett. 389 39-42 

[105] Casida M E 2005 /. Chem. Phys. 122 0541 1 1-9 

[106] Tozer D J and Handy N C 2000 Phys. Chem. Chem. Phys. 2 2117-21 

[107] Hirata S and Head-Gordon M 1999 Chem. Phys. Lett. 302 375-82 

[108] Tatarczyk K, Schindlmayr A and Scheffler M 2001 Phys. Rev. B 63 235106 

[109] Krieger J B, Li Y and Iafrate G J 1992 Phys. Rev. A 45 101-26 

[110] Kim Y and Gorling A 2002 Phys. Rev. B 66 35 1 14 

[111] Lein M, Gross E K U and Perdew J P 2000 Phys. Rev. B 61 13431-7 

[112] Burke K, Petersilka M and Gross E K U 2002 Recent Advances in Density Functional Methods 3 67-79 

[113] Corradini M, Del Sole R, Onida G and Palummo M 1998 Phys. Rev. B 57 14569-71 

[1 14] Moroni S, Ceperley D M and Senatore G 1995 Phys. Rev. Lett. 75 689-92 

[115] Richardson C F and Ashcroft N W 1994 Phys. Rev. B 50 8170-81 

[116] Hubbard J 1958 Proc. R. Soc. Lond. 243 336-52 

[117] Vashishta P and Singwi K S 1972 Phys. Rev. B 6 875-87 

[118] Vashishta P and Singwi K S 1972 Phys. Rev. B 6 4883(E) 

[119] Utsumi K and Ichimaru S 1980 Phys. Rev. B 22 5203-12 

[120] Dobson J F 1994 Phys. Rev. Lett. 73 2244-7 

[121] Dreizler R M and Gross E K U 1990 Density Functional Theory (Berlin: Springer) 

[122] Lipparini E, Stringari S and Takayanagi K 1994 /. Phys.: Condens. Matter 6 2025-30 

[123] Mermin N D 1970 Phys. Rev. B 1 2362-3 

[124] Kohn W 1996 Phys. Rev. Lett. 76 3168-71 

[125] Palummo M, Onida G, Del Sole R, Corradini M and Reining L 1999 Phys. Rev. B 60 1 1329-35 

[126] Gunnarsson O, Jonson M and Lundqvist B I 1979 Phys. Rev. B 20 3136-64 

[127] Olevano V, Palummo M, Onida G and Del Sole R 1999 Phys. Rev. B 60 14224-33 

[128] van Leeuwen R 1996 Phys. Rev. Lett. 76 3610-3 

[129] Hedin L 1965 Phys. Rev. 139 A796-823 

[130] Strinati G 1988 Riv. Nuovo Cimento 11 1-86 

[131] Sham L J and Schluter M 1983 Phys. Rev. Lett. 51 1888-91 

[132] Bruneval F, Sottile F, Olevano V and Reining L 2006 /. Chem. Phys. 124 1441 13 

[133] Kim Y and Gorling A 2002 Phys. Rev. Lett. 89 096402 

[134] Sharp R T and Horton G K 1953 Phys. Rev. 90 317 

[135] Talman J D and Shadwick W F 1976 Phys. Rev. A 14 36-40 

[136] Gorling A 1996 Phys. Rev. B 56 7024-9 

[137] Gorling A 1999 Phys. Rev. B 59 10370(E) 

[138] Stadele M, Moukara M, Majewski J A, Vogl P and Gorling A 1999 Phys. Rev. B 59 10031-43 

[139] Kotani T 1994 Phys. Rev. B 50 14816-21 

[140] Kotani T 1995 Phys. Rev. B 51 13903-4(E) 

[141] Kotani T 1995 Phys. Rev. Lett. 74 2989-92 

[142] Kotani T and Akai H 1996 Phys. Rev. B 54 16502-14 



Time-dependent density-functional theory for extended systems 



407 



[143] Stadele M, Majewski J A, Vogl P and Gorling A 1997 Phys. Rev. Lett. 79 2089-92 

[144] Magyar R J, Fleszar A and Gross E K U 2004 Phys. Rev. B 69 0451 1 1 

[145] Perdew J and Levy M 1983 Phys. Rev. Lett. 51 1884-7 

[146] Tokatly I V and Pankratov O 2001 Phys. Rev. Lett. 86 2078-81 

[147] Brosens F, Devrese J T and Lemmens L F 1980 Phys. Rev. B 21 1363-79 

[148] Aulbur W G, Stadele M and Gorling A 2000 Phys. Rev. B 62 7121-32 

[149] Griming M, Marini A and Rubio A 2006 /. Chem. Phys. 124 154108 

[150] Bruneval F, Sottile F, Olevano V, Del Sole R and Reining L 2005 Phys. Rev. Lett. 94 186402 

[151] Del Sole R, Reining L and Godby R W 1994 Phys. Rev. B 49 8024-8 

[152] Streitenberger P 1984 Phys. Status Solidi b 125 681-92 

[153] Streitenberger P 1984 Phys. Lett. A 106 57-60 

[154] Del Sole R, Adragna G, Olevano V and Reining L 2003 Phys. Rev. B 67 045207 

[155] Casida M E 1995 Recent Advances in Density Functional Methods vol 1 ed D P Chong (Singapore: World 
Scientific) p 155 

[156] Lautenschlager P, Garriga M, Vina L and Cardona M 1987 Phys. Rev. B 36 4821-30 

[157] Tokatly I V, Stubner R and Pankratov O 2002 Phys. Rev. B 65 1 13107 

[158] Hedin L and Lundqvist S 1969 Solid State Physics vol 23 (New York: Academic) pp 1-181 

[159] de Boeij P L, Kootstra F, Berger J A, van Leeuwen R and Snijders G 2001 J. Chem. Phys. 115 1995-9 

[160] van Faassen M, de Boeij P L, van Leeuwen R, Berger J A and Snijders J G 2003 /. Chem. Phys. 118 1044-53 

[161] Marini A, Del Sole R and Rubio A 2006 Time-Dependent Density Functional Theory (Lecture Notes in Physics 
vol 706) ed M A L Marques et al (Berlin: Springer) pp 301-16 (chapter 20) 
For the case of surfaces, see Pulci O, Marini A, Palummo M and Del Sole R 2007 Phys. Rev. B submitted 

[162] Varsano D, Marini A and Rubio A 2006 submitted 

[163] Maitra N T, Sousa I and Burke K 2003 Phys. Rev. B 68 045109 

[164] Dhara A K and Ghosh S K 1987 Phys. Rev. A 35 442-4 

[165] D'Agosta R and Vignale G 2005 Phys. Rev. B 71 245103 

[166] Vignale G 2004 Phys. Rev. B 70 201 102(R) 

[167] Vignale G and Kohn W 1996 Phys. Rev. Lett. 77 2037-40 

[168] Vignale G, Ullrich C A and Conti S 1997 Phys. Rev. Lett. 79 4878-81 

[169] Vignale G 2001 Int. J. Mod. Phys. 15 1714-23 

[170] Ullrich C A and Vignale G 2001 Phys. Rev. Lett. 87 037402 

[171] Ullrich C A and Vignale G 2002 Phys. Rev. B 65 245102 

[172] Ullrich C A and Vignale G 2004 Phys. Rev. B 70 239903(E) 

[173] van Faassen M, de Boeij P L, van Leeuwen R, Berger J A and Snijders J G 2002 Phys. Rev. Lett. 88 186401 

[174] Ghosez P, Gonze X and Godby R W 1997 Phys. Rev. B 56 12811-7 

[175] Sottile F, Karlsson K, Reining L and Aryasetiawan F 2003 Phys. Rev. B 68 205112 

[176] Godby R W and Sham L J 1994 Phys. Rev. B 49 1849-57 

[177] Aulbur W G, Jonsson L and Wilkins J W 1996 Phys. Rev. B 54 8540-50 

[178] Vignale G and Kohn W 1998 Electronic Density Functional Theory: Recent Progress and New Directions 

ed J Dobson et al (New York: Plenum) p 199 

[179] Benedict L X, Shirley E L and Bohn R B 1998 Phys. Rev. B 57 R9385-7 

[180] Rowe J and Aspnes D 1970 Phys. Rev. Lett. 25 162-5 

[181] Rowe J and Aspnes D 1970 Phys. Rev. Lett. 25 979 (E) 

[182] Martin R M, van Vechten J A, Rowe J E and Aspnes D E 1972 Phys. Rev. B 6 2500-2 

[183] Schwabe N F and Elliott R J 1996 Phys. Rev. B 53 5318-29 

[184] Egri I 1985 Phys. Rep. 119 363-402 and references therein 



